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Abstract 


Extensions of the Hartree-Fock-Bogoliubov theory are worked out 
which are tailored for, (i) the consistent evaluation of fluctuations 
and correlations and (ii) the restoration through projection of broken 
symmetries. For both purposes we rely on a single variational princi¬ 
ple which optimizes the characteristic function. The Bloch equation 
is used as a constraint to define the equilibrium state, and the trial 
quantities are a density operator and a Lagrangian multiplier matrix 
which is akin to an observable. The conditions of stationarity are re¬ 
spectively a Schrodinger-like equation and a Heisenberg-like equation 
with an imaginary time running backwards. General conditions for 
the trial spaces are stated that warrant the preservation of thermody¬ 
namic relations. The connection with the standard minimum principle 
for thermodynamic functions is discussed. When the trial spaces are 
chosen to be of the independent-quasi-particle type, the ensuing cou¬ 
pled equations provide an extension of the Hartree-Fock-Bogoliubov 
approximation, which optimizes the characteristic function. An ex¬ 
pansion of the latter in powers of its sources yields for the fluctuations 
and correlations compact formulae in which the RPA kernel emerges 
variationally. Variational expressions for thermodynamic quantities or 
characteristic functions are also obtained with projected trial states, 
whether an invariance symmetry is broken or not. In particular, the 
projection on even or odd particle number is worked out for a pairing 
Hamiltonian, which leads to new equations replacing the BCS ones. 
Qualitative differences between even and odd systems, depending on 
the temperature T, the level density and the strength of the pairing 
force, are investigated analytically and numerically. When the single¬ 
particle level spacing is small compared to the BCS gap A at zero 
temperature, pairing correlations are effective, for both even and odd 
projected cases, at all temperatures below the BCS critical tempera¬ 
ture Tc. There exists a crossover temperature Tx such that odd-even 
effects disappear for Tx < T < Tc. Below Tx, the free-energy differ¬ 
ence between odd and even systems decreases quasi-linearly with T. 
The low temperature entropy for odd systems has the Sakur-Tetrode 
form. When the level spacing is comparable with A, pairing in odd 
systems takes place only between two critical temperatures, exhibiting 
thus a reentrance effect. 
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1 Introduction 


Let us consider a finite physical system in thermodynamic equilibrium. We 
want to evaluate variationally the thermodynamic functions of this system 
along with the expectation values, fluctuations and correlations of some set 
of observables Q^. The equilibrium state of the system is dehned from our 
knowledge of conserved quantities such as the number of particles, the energy, 
the total angular momentum, the parity. However, because our system of 
interest is hnite, these equilibrium data should be handled with care. Indeed, 
depending on the circumstances, their values can be given in two different 
ways, either with certainty or statistically, and this allows us to classify 
the conserved quantities into two types. Take for instance a hnite, closed 
system thermalized with a heat bath. Among the variables that determine 
its state, the particle number N is given with certainty ; we shall call N a 
conserved quantity of type 1. On the other hand, the energy is known only 
statistically as only its expectation value is hxed; let us call it a conserved 
quantity of type 11. The Gibbs argument, which identihes a heat bath with 
a large set of copies of the system, entails that this system is in a canonical 
Boltzmann-Gibbs equilibrium characterized by a Lagrangian multiplier (3 
dehning the inverse temperature. Apart from N, our information is specihed 
by (3 or equivalently by the expectation value of the Hamiltonian. Take now 
an open system which communicates with heat and particle baths. In this 
case, both the energy and the particle number are conserved quantities of 
type 11, and the Gibbs argument leads to the grand canonical distribution. 
For other equilibrium data one has likewise to recognize, according to the 
circumstances, whether they are of type 1 or type 11. Note that, aside from 
the truly conserved quantities, the type 11 equilibrium data may also include 
nearly conserved quantities (or order parameters) which are here assumed to 
have all been identihed beforehand. 

For macroscopic systems, this distinction between types 1 and 11 is imma¬ 
terial because the relative fluctuations that occur for observables of type 11 
are negligible in the thermodynamic limit. Ganonical and grand canonical 
distributions then become equivalent. When dealing with finite systems, 
however, it is important to specify carefully which statistical ensemble is 
relevant to the physical situation under consideration. 

This problem is encountered in various systems of current interest. In nu¬ 
clear physics, depending on the experimental conditions and/or the accuracy 
of the description, we may have to describe nuclei having either a well-dehned 
or a fluctuating particle number, angular momentum or energy. Similar situ¬ 
ations may occur for atoms or molecules. Another example is that of metallic 
or atomic clusters which can be prepared with a given particle number []^. 
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Let us also mention mesoscopic rings for which a proper description of the 
conductivity properties requires that there be no fluctuations (i.e., a canon¬ 
ical distribution) in the number of electrons|0]. Likewise, the description of 
mesoscopic superconducting islands 0 or metallic grainsdemands a distri¬ 
bution with a given parity of the number of electrons. The recent discovery of 
Bose-Einstein condensates in dilute alkali atomic gases is currently arousing 
great interest P| ; in this case also, the theoretical formulation requires us to 
use states with well-defined numbers of particles. 

A standard formal method to assign a density operator Z) to a system 
characterized by data of both types proceeds as follows^. The type I data, 
referring to some conserved observables, determine a Hilbert space, or a 
subspace Tii (for instance the subspace with a given number of particles) in 
which the density matrix D operates. On the other hand, type II data are 
the expectation values (0^) of the remaining conserved observables which 
act in Til. As a consequence, they impose on D the constraints 

(Ofc) = Tn {ikD/Tii D , (1.1) 


where the trace Tri is meant to be taken on the Hilbert subspace TZi. As 
we shall see, it is convenient to leave D unnormalized or, equivalently, not 
to include the unit operator among the set 0^. Within the subspace 7Zi the 
state D is determined, up to a normalization factor, by maximizing the von 
Neumann entropy 


S = 


Tri DhiD 
Tri b 


+ liiTri b 


( 1 . 2 ) 


under the constraints (|1.1|) . Associating a Lagrangian multiplier with each 
datum (Hfc), one recovers the generalized Boltzmann-Gibbs distribution 


b (X e 


(1.3) 


where and (Clk) are related to each other through 


i^k) — 


d 


dXi 


InTi’i e 


(1.4) 


In what follows, we take the form ( p..3|) of the density operator for granted, 
irrespective of the mechanism that led to it. Equilibrium distributions of the 
form ( |1.3|) can be the outcome of some transport equation, the detailed nature 
of which is beyond the scope of this work. The form (|1.3|) can also be justi¬ 
fied by various theoretical arguments. One of them is the above-mentioned 
maximum entropy criterion which has its roots in information theory 0. Al¬ 
ternatively, the Gibbs argument provides a direct statistical derivation of the 
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Boltzmann-Gibbs distribution. An extension to non-commuting observables 
Clk, relying as in the classical case on Laplace’s principle of indifference, has 
been worked out in ref. [Q. 

Formally, the assignment of the density operator (|1.3|) to the system solves 
our problem since the thermodynamic functions, expectations values, fluctua¬ 
tions, correlations of the observables and all their cumulants are conveniently 
obtained from the characteristic function 

(p(0 = hiTn i(0^ , (1.5) 

where the operator A(^), depending on the sources associated with the 
observables Q 7 , is dehned as 

i (0 = exp(-^^^Q^) . ( 1 . 6 ) 

7 

Indeed, by expanding ([T75|) in powers of the sources, 

V>(?) = »’(0) - E + 5 E ■ (1.7) 

7 7(5 

one gets at hrst order the expectation values (Q 7 ) of the observables Q^. At 
second order appear their correlations 

= \{Q'yQ5 + QsQ'y) — {Q'ftiQs) (1.8) 

and fluctuations = (^ 77 , and in higher orders appear the cumulants of 
higher rank. Moreover, since we chose not to enforce the normalization of D, 
the generalized thermodynamic potential is given by the zeroth-order term : 

(p(0) = InTri b = InTri e" = S-Y, Afc(Gfc) . (1.9) 

k 

Unfortunately, in most physical cases, the calculation above cannot be 
worked out explicitly. Indeed, the statistical operator (|1.3|) is tractable only if 
all the operators VLk are of the single-particle type. For systems of interacting 
particles, this is no longer the case since, in general, the data (Gfc) include the 
average energy {H). Recourse to some approximation scheme then becomes 
unavoidable. Many approaches start with the replacement of the density 
operator ( [1.3|) by one with the form 

Poce"* , ( 1 . 10 ) 

where M is a single-particle operator J2ij MijO^iOj, or more generally a single 
quasi-particle operator which includes also pairs and 0 * 0 ^. (In this 
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article, we shall mainly deal with systems of fermions ; for the boson case, 
terms linear in the operators a'*'* and Oj should be added to M.) 

However, the form ( p..lU| ) entails two difficulties, both of which will be 
adressed in this article. 

(i) The hrst difficulty concerns the optimum selection of the independent- 
quasi-particle approximate state (JL 13). One often uses a variational criterion, 
and this is an approach that we shall follow also. [Conventional many-body 
perturbation theory is not suited to the problems we have in mind since no 
small parameter is available.] The choice of this criterion, however, raises 
a crucial question. A standard option, which leads to the Hartree-Fock 
(or Hartree-Fock-Bogoliubov) solution, consists in choosing the minimiza¬ 
tion of the free energy or, at zero temperature, of the energy. This procedure 
is clearly adapted to the evaluation of the thermodynamic quantities that 
are given by </9(0), the value of Eq. (|^) when the sources vanish (see 
Eq. (|1.9|)). However, the resulting approximate state is not necessarily suited 
to an optimization of the characteristic function (p(^) for non-zero values of 
the sources For restricted trial spaces such as (|1.10|) (or (|1.11|) ), a vari¬ 
ational method designed to evaluate (p(0) rather than (p(0), is expected to 
yield an optimal independent-particle density operator which, owing to its 
dependence on the intensity of the sources, is better adapted to our purpose 
than the (^-independent) Hartree-Fock solution. 

(ii) The second difficulty is more familiar. Let us for instance consider 

a system for which the particle number is of type I and the energy of type 
H. It should be described by a canonical equilibrium distribution D oc 
built up in the A^-particle Hilbert space Tii. Since our system is hnite, 
the results are expected to differ from those of a grand canonical equilib¬ 
rium. However, perturbative and variational approaches are often worked 
out more conveniently in the entire Fock space Ti. In particular, even an 
independent-particle state such as (|1.10|) where M = 'ffij resides 

in this Fock space 7Y. Thus the simple Ansatz (|1.1U|) requires leaving the 
space TLi to which the exact canonical state belongs; as a consequence, the 
standard variational treatment introduces spurious components into the trial 
state. Unphysical thermodynamic fluctuations of the particle-number thus 
arise from the form of the trial state (|1.1CI|) , which would be better suited 
to a grand canonical than to the canonical equilibrium being approximated. 
These fluctuations should be eliminated, and it is natural to improve the 
approximation by introducing a projection onto states with the right number 
of particles. This is achieved by replacing Eq. (|1.1CI|) by the more elaborate 
Ansatz 

V(xPn e-^ Pn , 


(1.11) 





















where Pn is the projection onto the subspace Hi with the exact particle 
number N. The characteristic function ([1.5|) can now be evaluated by taking 
the trace Tr in the full Fock space H. Likewise, when the exact state has 
a non-fluctuating angular momentum or spatial parity, any trial state of 
the form ( LlOj ) introduces thermodynamic fluctuations for these quantities, 
and adequate projections P are required so as to suppress these fluctuations. 
However, if the operator M commutes with P, it is not necessary in Eq. (|1.11|) 
to insert the projection P on both sides of the exponential. 

Projections turn out to be very useful in another circumstance, namely 
when a symmetry is broken by the operator M. Since in this case M does 
not commute with P, the projection should appear on both sides as in (|1 .1 1|) . 
Although projection matters mostly for finite systems, let us hrst recall, for 
the sake of comparison, the situation in the inhnite systems. 

Consider a set of interacting particles which, in the thermodynamic limit, 
presents a spontaneously broken symmetry. One can think of a Bose liquid 
or a superconducting electron gas where the A^-invariance is broken. In such 
a system an essential property of the exact statistical operator D is the 
long-range order displayed by the associated n-point functions. For an in¬ 
teracting Bose fluid at low temperature, if we denote by V^^(r) the creation 
operator at the point r, expectation values over D such as ('^'*'(ri)-^(r 2 )), 
(V’^(ri)'^'*'(r 2 )'^(r 3 )' 0 (r 4 )) factorize in terms of a single function (p(r) as <p(ri)*<p(r 2 ), 
<p(ri)*<p(r 2 )*</ 7 (r 3 )<p(r 4 ) when the various points all move apart from one an¬ 
other. On the other hand, since for any size of the system its exact density 
operator D commutes with N, as does the Hamiltonian H, the quantity 
('0(r)) vanishes. In order to give a meaning to the factors </ 3 (r), a stan¬ 
dard procedure consists of adding to H small sources which do not commute 
with N: for instance terms in a* and a^i for bosons, or in a^ia^j and OiOj 
for fermions. The resulting equilibrium density operator Db displays there¬ 
fore an explicit broken invariance. Despite the fact that such sources do 
not exist in nature, the remarkable properties of Db make it an almost in¬ 
evitable substitute for D. Indeed, in the limit when one hrst lets the size 
of the system grow to inhnity and then lets the sources tend to zero, Db is 
equivalent to D as regards all the observables that commute with iV, such 
as V’^(ri)^/’(r 2 ) or ' 0 l(ri)'^l(r 2 )'^(r 3 )'^(r 4 ); however, the expectation value 
(V’(r))B over Db does not any longer vanish and is interpreted as an order 
parameter. Furthermore, Db satishes the clustering property for all opera¬ 
tors, even those which do not commute with iV; in particular, we now have 
(^/’■f(ri)' 0 (r 2 ))B ^ (' 0 Hri))B(' 0 (r 2 ))B when |ri - r 2 | ^ oo^so that ('0'^(r))B 
can be identihed with <p(r). The clustering property of Db as well as the 
equivalence of D and Db for the physical observables were recognized for the 
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pairing of electrons as early as the inception of the BCS theory^. Thus, 
for sufficiently large systems, the same physical results are obtained not only 
for different canonical ensembles, but also with or without the inclusion of 
sources that break invariances. 

The situation is different for a hnite system governed by the same inter¬ 
actions as in the infinite case above, under similar conditions of density and 
temperature. Here again the expectation value ('0(r)) is zero in the exact 
equilibrium state, but now no analogue to D-q can be devised. Indeed, we 
can no longer define an order parameter by means of factorization of n-point 
functions since the relative coordinates cannot grow to infinity. Nor can 
we introduce sources which would produce non-vanishing expectation values 
('0(r))B without changing the physical n-point functions, as this requires the 
sources to tend to zero only after the thermodynamic limit has been taken. 
Even the translational invariance of a nucleus, an atom or a cluster cannot 
be broken in the exact ground state or in thermal equilibrium states if the 
Hamiltonian includes no external forces and commutes with the total mo¬ 
mentum. Returning to the iV-invariance, for a nucleus with pairing or for a 
superconducting mesoscopic island or metallic grain, it is no longer without 
physical consequences to break the e**^^ gauge invariance. This not only in¬ 
troduces spurious off-diagonal elements in N for the density matrix Z), but 
also affects observable matrix elements in the diagonal blocks in N. 

Thus, the breaking of invariances does not have the same status for hnite 
and inhnite systems. In the latter case, it stands as a conceptual ingredi¬ 
ent of the theoretical description which correctly implements the clustering 
property of the n-point functions. For hnite systems it is also essential, but 
now as a basic tool for building approximation schemes which economically 
account for some main features of the correlations. Both in perturbative 
and variational treatments, the starting point is most often an approximate 
density operator of the form ( p..lO| ) which does not commute with conserved 
observables such as N. This yields tractable and sometimes remarkably ac¬ 
curate approximations, for which the clustering of the correlation functions 
arises naturally from Wick’s theorem. Remember the success of the BCS 
theory in describing pairing correlations betweens nucleons in hnite nuclei. 
(One may wonder why nuclear physicists, who had acknowledged the exis¬ 
tence of pairing ehects long before 1957, did not anticipate the BCS theory ; 
it may be that, too obsessed by the conservation of the particle number, they 
were inhibited in developping models which would break the iV-invariance.) 

Nonetheless, more realistic descriptions of hnite systems require the restora¬ 
tion of the broken invariances. A state of the form ( 1.1U| ) involves huctuations 
in the particle number that must be eliminated, whether they arise from a 
grand canonical description, or from unphysical oh-diagonal elements of N 
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connecting Hilbert spaces with different particle numbers. A natural pro¬ 
cedure to account for pairing correlations and at the same time enforce the 
constraint on the particle number consists in using a projection. The use of 
Eq. (|1 .1 1|) arises naturally in this context. Likewise, breaking translational or 
rotational invariance introduces spurious fluctuations of the linear or angular 
momentum, as well as spurious off-diagonal elements in the trial density op¬ 
erator (ll.lOP . To remedy these defects we need again a trial density operator 
of the form ( p..llj ) involving the appropriate projection P on both sides. 

The discussion above pertains to the breaking of a conserved quantity of 
type I, that is, one which is known with certainty , such as in a canonical 
ensemble. For a quantity of type II, known only on average, such as N in 
a grand canonical ensemble, the invariance breaking again introduces off- 
diagonal elements of N, which are spurious for a finite system. One can 
suppress them, while letting N free to fluctuate, by means of a different type 
of projection, exhibited below in Eq. (|5.11|) . 

The problem of the restoration of broken symmetries has already had a 
long history. Within the nuclear context, in the zero-temperature limit, a 
general review of the litterature (up to 1980) may be found in |]^. For a 


recent work about projection at finite temperature and its connection with 
the static path approximation see Ref. [^, where more references can be 


found. An earlier reference on parity projection is 


For a Monte-Carlo 


approach to the temperature dependence of pair correlations in nuclei, see 
Ref.[|^]. In condensed matter physics, recent experimental studies about 
isolated mesoscopic metallic islands 0 or small metallic grains ^ have moti¬ 
vated several theoretical works about modified BCS theories with well defined 
number parity||I^-pm. 

The two problems (i) and (ii) that we have stated above will be worked 
out within a general variational setting described in Sect. 2. We shall write an 
action-like functional, Eq. ( p.4|) , which yields as its stationary value the char¬ 
acteristic function (|1.5|). This corresponds to the static limit of a variational 


principle proposed elsewhere to evaluate multi-time correlation functions 

123, 


The approach uses a method 


which allows the construction of vari¬ 


ational principles (and variational approximations) tailored to the optimal 
evaluation of the quantity of interest (here the characteristic function). This 
systematic and very general method encompasses several known variational 
principles, including the Lippmann-Schwinger formalism (the quantity 
of interest being the S'-matrix). An important ingredient is the incorporation 
in the trial functional of the Lagrangian multipliers associated with the pri¬ 
mary equations, in our case the (symmetrized) Bloch equation determining 
the canonical density operator ( |1.3|) . As a result, the number of variational 
degrees of freedom is doubled. Here, the variational quantities entering our 
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functional (|2.4|) consists of two matrices : one akin to a density operator and 
the other to the observable dehned in (|1.6|) . Within this framework, one 
can state some general conditions on the trial spaces that ensnre the preser¬ 
vation of the thermodynamic identities by the variational approximations. 

Starting in Sect. 3, we apply the variational principle of Sect. 2 to sys¬ 
tems of fermions for which pairing effects are significant. Sect. 3 (as well 
as Sect. 4) deals with grand canonical eqnilibrinm. The two variational op¬ 
erators introdnced in Sect. 2 are taken as exponentials ( p.4p of single-qnasi- 
particle operators and the resnlting expression of the approximate action-like 
fnnctional is given in Sect. 3.2. The equations which express the stationarity 
and determine the optimal characteristic fnnction, are derived in Sect 3.3. 
They couple, with mixed boundary conditions, the variational parameters 
characterizing the two trial matrices and they generalize the Hartree-Fock- 
Bogolinbov (HFB) approximation!^ . 

In Sect. 4, the coupled equations obtained in Sect. 3 are expanded in 
powers of the sonrces ^.y. The HFB approximation is recovered at the zeroth 
order, while the first order gives back the usual HFB formula for the average 
values of observables (Sect. 4.1). Fluctuations and correlations of observables 
are provided by the second order. In the ensning formulas for conserved and 
non-conserved observables, the RPA (random-phase approximation) kernel 
plays a central role. In Sect. 4.3.3, throngh a diagrammatic analysis, we 
show the correspondence between onr resnlts and the snmmation of the RPA 
(bnbble) diagrams. 

Sect. 5.1 addresses the projection problem stated in (ii) above. The fnnc¬ 
tional and the conpled eqnations that result from variational Ansatze of the 
type (|1.11|) are derived in Sects. 5.2 and 5.3, respectively. Special attention 
is devoted to the case where the variational approximation does not break 
the P-invariance (Sect. 5.4). 

Sect. 6 specializes the formalism of Sect. 5 to the projection on even or 
odd particle nnmber. The more detailed stndy of this relatively simple case 
is motivated by its relevance to the description of heavy nnclei and to the 
interpretation of recent experiments on isolated supercondncting islands]^ 
or small metallic grains]^. Equations are obtained in Sect. 6.3 which should 
replace the BCS ones. Ensning results are compared to those of Ref. . In 
Section 6.4 we analyze some limiting situations with a special emphasis on 
the low temperatnre limit. Numerical applications are presented in Sect. 6.5; 
they illustrate physical situations enconntered in nnclei, metallic islands or 
grains. 

In the last Section we snmmarize onr main results and present some 
snggestions for other applications and for extensions. 

Appendix A introdnces onr notation for the Lionville representation of 
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the reduced fermion quasi-particle space. In particular, it gives a proof of the 
factorization of the RPA kernel (used in Sect. 4) into the stability matrix and 
the matrix whose elements are the constants associated with the Lie-Poisson 
structure of the HFB equations. 

In Appendix B, using the Liouville-space notation, we reexpress the equa¬ 
tions derived in Sect. 5 in a fashion which exhibits their formal similarity with 
those obtained in Sect. 3. 

The length of this paper is partly a consequence of our wish to allow sev¬ 
eral reading options. The contents of Sects. 5 and 6 are largely independent 
of that of Sect. 4. The reader primarily interested by the variational pro¬ 
jection method can therefore bypass Sect.4, except for the brief Sect. 4.1. If 
he is mostly concerned by the projection on even or odd particle-number, 
and ready to admit a few formulas demonstrated in Sect. 5, he can di¬ 
rectly jump to Sect. 6. If his interest is focused upon the improvements 
over BCS introduced by parity-number projection, he may even begin with 
Sect. 6.3. Alternatively, the reader mainly interested by the variational eval¬ 
uation of correlations and fluctuations in the grand canonical formalism can 
skip Sects. 5 and 6. 
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2 The Variational Principle 

Our purpose in this Section is to write a variational expression adapted to 
the evaluation of the characteristic function (|1.5|) , namely 

exp^(0=Trii(0l) . (2.1) 

The operator A, dehned in Eq. o, involves the observables of interest Q-y 
and depends on the associated sources The density operator D describes 
thermodynamic equilibrium. It is an exponential of the conserved quantities 
Clk of type II that are only given on average (Eq. (|1.3|) ). From now on we 
assume that the Hamiltonian H is one of these conserved quantities, with (3 
as its conjugate variable. We thus write D in the form 

D = e-f^^ , l3k = Y.\khk . ( 2 . 2 ) 

k 

We have not introduced the unit operator among the set {12^.}. Then </9(0) 
does not vanish but is identihed (Eq. (|1.9|) ) with the logarithm of the partition 
function Tri , for instance of the canonical partition function for K = H 
or of the grand canonical partition function for K = H — fiN. Thus, our 
variational principle will yield not only the expectation values and cumulants 
of the set {Qy}, but also the thermodynamic quantities. As regards the 
conserved quantities of type I, they dehne the Hilbert space Hi over which 
the trace Tri is taken in Eq. ( |2.1|) . According to the circumstances. Hi can 
be the Fock space, as in Sects. 3 and 4, or a more restricted space. Sects. 5 
and 6 resort in principle to the latter case; nonetheless, the introduction, as 
in Eq. ( p..ll|) , of a projection onto Hi will allow us to use again the trace Tr 
in the full Fock space. 


2.1 The Action-like Functional and the Stationarity 
Conditions 


We wish to replace the density operator (|2.2| ) by a variational approximation 
P(/3), since D itself is untractable. With this in mind, we note that D = e“^^ 
can be constructed as the solution of the Bloch equation, or rather of 
its symmetrized form 


dT) 1 - - - -, 

— + -[KV + VK] =0 
du 2 


(2.3) 


where T’(m) is a u-dependent operator {0 < u < (3) subject to the initial 
condition I?(0) = 1. In accordance with the general method |^, ^ devised 
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to evaluate variationally a quantity of interest, here ( p.lD , we regard the 
equation (|2.3|) as a set of constraints which impose the canonical form V{f3) = 
, and we associate with it a Lagrangian multiplier A{u) which is also a u- 
dependent operator. The characteristic function O is then the stationary 
value of the action-like functional 


I{V{u), A{u)} = Tri AV{I3) 

- du Tri A{u) ( + hkV{u) +V{u)k] 

Jo \ du 


(2.4) 


under arbitrary variations of the trial operators A{u) and P(m), the latter 
being subject to the boundary condition 


V{0) = i 


(2.5) 


The specihcs of our problem enter the functional (|2.4|) through the tem¬ 
perature 1//9, the operator k and through the operator A (a functional of 
the observables). 

A variational approximation for exp93(,^) will be generated by restrict¬ 
ing the trial class for A{u) and T’(m) and by looking for the corresponding 
stationary value of (|2.4|) . The error will be of second order in the difference 
between the optimum trial operators and the exact ones. An essential fea¬ 
ture of the method is the doubling of the number of variational degrees of 
freedom, which include not only the trial state T){u) but also the Lagrangian 
multiplier operator A{u). This doubling is beneficial only when approxima¬ 
tions are made; indeed, for unrestricted variations of A{u) the stationarity 
conditions on A{u) alone reduce to the Bloch equation ( |2.3|) and are sufficient 
to ensure that T{V^ A} yields the exact value for exp93(,^). 

When V{u) and A{u) are restricted to vary within some subset, the sta¬ 
tionarity conditions of the functional X{I?, A} with respect to A{u) are read 


directly from Eq. (|2.4|) : 


Tri SA 



+ -[kV + VK] 


0 


( 2 . 6 ) 


As regards the variation with respect to T>{u), an integration by parts allows 
us to rewrite the functional (^.41) in the form 


I{V{u), A{u)} = Tri ^(0) + Tri X(/?)(i - A{(3)) 


d- du Tri Xfn) 
Jo 


' dA{u) 
du 


-[A{u)k + kA{u)] 


(2.7) 
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From Eq. (E3). for 0 < M < /3, one obtains the stationarity conditions: 


Tr, iS® - \[Ai< + = 0 , 


( 2 . 8 ) 


and for n = /?: 

i:ii5V{l3){A-A{l3)) = Q . (2.9) 

In Eqs. ( |2.6| ) and ( p.8|) we have omitted the w-dependence of "D, of A, and 
of the allowed variations hP and 6A. 

By construction of the functional (|2.4|) , Eq. (p.6| ) gives back the sym¬ 
metrized Bloch equation (|2.3|) for unrestricted variations of A. For unre¬ 
stricted variations of P, Eqs. ( p.8|) and (|^) yield (in the Hilbert space Tii) 
the Bloch-like equation for A{u) 


d i 1 - - - - 

'^--[AK + KA]=Q 
du 2 

in the interval 0 < u < f3, and the boundary condition 

A{p) = A 


( 2 . 10 ) 


( 2 . 11 ) 


at the hnal value u = 13. The variable u is thus evolving backward from (3 
towards 0 in Eq. (|2.10|) , the solution of which is 


A{u) = exp[i(M — /d)iF]y4exp[i(M — I3)K] 


( 2 . 12 ) 


The exact equation ( p.l0| ), which is analagous to a backward Heisenberg 
equation in the “time” n, merely duplicates the Schrodinger-like equation 
(|2.8|). However, in any physical application, the trial spaces for T>{u) and 
A{u) have to be restricted to make the calculations tractable. Then, in 
general, Eqs. ( [^.6| ) and ( [^.8| ) are coupled and must be solved simultaneously. 
Moreover, the boundary conditions (|2.5| ) on T>{d) and (|2.9|) on A{(3) imply 
that Eqs. 


and (|2.8|) should be solved in the forward and backward 
directions in u, respectively. 


2.2 General Properties 

The stationary value Xg of the functional ( 12.4|) depends on an ensemble of 
parameters r]. Among them are obviously the source parameters [the 
dependence is through the boundary condition (|2.9|) ] and the temperature 
1//3. Other parameters a; may occur in the operator K, as when it includes 
an external field —ujJ. In particular, for a grand canonical ensemble, where 
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K = H — the chemical potential /x is regarded as one of the parameters 
ct). The symbol rj will therefore stand for the set [3, ut}. 

We shall take advantage of a general property of variational principles 
that we now state in onr notation. The stationary valne Ts{ri) which is 
songht, 

Uv) ^ Xs{i(0, /9, K} = 1{K K v) , (2.13) 

depends on the parameters r] both directly and throngh the optimal soln- 
tion determined within onr variational class by the stationarity 

conditions 5X{V, A, = 0 ^ind 5X{'D, A, t]}/6A = 0. Owing to these 

conditions, the fnll derivative of Xg with respect to rj rednces to its partial 
derivative evalnated at the stationary point : 



dr] 


■X{V, A, T]} 


'D—'D^ 

A=Ar] 


(2.14) 


This relation holds not only for the exact solntion bnt also within any trial 
class for A and V. 

For any variational space, the derivatives dT>/du and dA/du belong to 
the set of allowed variations 5V and 6A , respectively. Writing Eqs. ( ^TB| ) 
and ( p.8|) for these variations and snbtracting one eqnation from the other, 
we obtain along the stationary path the conservation law 

k [Ar^{u)Vn{u)+'br^{u)A^{u)] = 0 . (2.15) 


In the following, we shall always nse for A a trial space 
nnit operator and in which variations 5A proportional to 
allowed. We then End from Eq. (|2.6|) that the integrand 
(T^) vanishes, so that the stationary valne is given by 


that contains the 
A {6A oc A) are 
of the fnnctional 


XM = Tn AV,{(3) 


(2.16) 


as in the exact expression m- 

Moreover when, as will be the case in this work, the trial space of T> is 
chosen snch that variations 61) proportional to T) {61) oc X)) are allowed, 
another helpfnl property emerges (which of conrse holds also for the exact 
solntion). Indeed, snbstitnting 6A oc M in Eq. ( |2.(j|) and 6V oc X in Eq. (PTB|) 
and combining these eqnations, we obtain 

^Tri ArjikV^iu) = 0 . (2.17) 
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Using 5V (xV in Eq. we also find 

Tn AV^{(3) = Tn Ar,{(3)V^{(3) , (2.18) 

even when A does not belong to the variational space for A. The rela¬ 
tions (|2.17|) and (|2.18|) imply that the approximation ( p.l6|) for exp(p = 
Tri AVriif^) is equal to Tri for any value u, and in particular 

that 

g(p(0 = Tri A(0) . (2.19) 

Let us now use Eq. (|2.14|) when 7] stands for the source parameters ; 
we hnd : 


de 


K}=e 




( 2 . 20 ) 


This relation is valid for all values of the ^’s. According to Eq. o, in the 
limit where these ^’s vanish, Eq. (|2.20|) yields 


d(p 

d^ 


= iQ'y) = 


«=0 


Tn 

Xs{i, /3, k} 


( 2 . 21 ) 


where, as in the rest of this Section, the subscript .^ = 0 is shorthand for 
= 0, /3, a;}. Since it follows from Eq. (|2.16|) that 


one has 


J,{i, A = Tn X)^=o(/3) 

^ Tn Q^V^=o{{3) 


( 2 . 22 ) 


(2.23) 


Tri V^=o{P) 

This result tells us that the calculation of the expectation values of any op- 

]) for A = 1. Equation 


erator requires only the solution of Eqs. (|2.6|) and ([ 

( |2.23|) shows therefore that the density operator V^=q{/ 3) suited to the calcu¬ 
lation of the thermodynamic potentials is also suited to the variational eval¬ 
uation of the expectation values of the observables Q^. The knowledge of the 
^-dependence of T>^(/3) is only needed for the calculation of the higher-order 
cumulants. 

To hnd the derivative of Xg with respect to the inverse temperature (3, we 


rewrite the functional (|2.4|) as 


1 = 


duTTi\AV{u)6{u- (3) 


— A{u) 


'dViu 


(2.24) 


+-[KV{u)+V{u)K]\e{f3-u) 
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We then apply Eq. (|2.14|) with r] replaced by 13. Noting that d'D{u)/(3u\u=i3 
belongs to the variational class allowed for SV{j3) and using (|2.9|) , we obtain 


^ K IA,(I3)V,{P) + V,{P)APP)] 


(2.25) 


in which Arj{l3) can be replaced by A when the variational space for A con¬ 
tains A . When A = i, Eq. (|2.25|) becomes 


d/3 d/3 Ti'i®5.„(/3) ' ' 


(2.26) 


where we have used Eq. (|2.23| ) and introduced the notation E for the equi¬ 
librium energy in the canonical case {K = H), or the expectation value 
{H — fiN) in the grand-canonical case. Thus, in the framework of our vari¬ 
ational method ( for any trial space allowing 5A oc A), the derivative with 
respect to (3 of the approximate thermodynamic potential (A(0 = 1 ) pro¬ 


vides the same result as the expectation value of K obtained from Eq. (|2.23|) 
with Q = K (A(,^) = exp—We can also dehne an approximate ther¬ 
modynamics potential 

Fo = -/?-V(0) , (2.27) 


and an approximate entropy 


S = (3{E- Eg) = <,9(0) - 13 


d<,9(0) 

d(3 


(2.28) 


which, owing to Eq. 


, satisfy the relation 




OFg 
' dT 


(2.29) 


Our approximation is then consistent with the fundamental thermodynamic 
relations between the partition function, the temperature, the energy and 
the entropy. Note that S as dehned by Eq.(p.28|) is not necessarily related 
to V^=o{l3) by Eq. QT^) . 


If we now use the relation (|2.14|) where r] is identihed with the parameter 
uj of an external held —ojJ entering K, we have 


Tx, 

den dm 2. 


du Tri J [bn{u)AAu) + AAu)Vn{u)\ 


(2.30) 


When A = 1, i.e. when = 0, we shall see in Sect. 2.4 that under some 
rather general conditions (satished for instance in Sects. 4.1, 5.4.2, 6.2 and 


19 


















6.3), the products Vri{u)Ar,{u) = Ar,{u)Vr){u) do not depend on the variable 
u. Then Eq. (|2.30|) reduces to 


liL 

(3 do; 


InTn 


Tn JV^=o{l3) 
Tn P^=o(/9) 


(2.31) 


Here again, as for Eq. ( |2.26|) , our variational approximation is consistent 
with the thermodynamic identities expressing the expectation values of ob¬ 
servables J as derivatives of the thermodynamic potential. 

Let us hnally consider the case when the variational classes for T) and A 
are such that T>~^6V and 5AA~^ belong to a common set of operators. Then 
letting 5V = VX and dA = xA'm. Eqs. (^) and (|^) and subtracting one 
equation from the other, we hnd 


T:y,X['^+\[M>,k])=0 . (2.32) 

Likewise, if we let 5A = Ay and 51) = ¥1), we hnd 

Tril>(^^-i|Bi,A-]) = 0 . (2.33) 


For unrestricted variations, X and Y are arbitrary and Eqs. ( 2.32 - 0^) imply 
that AV and T>A both satisfy exact equations of the Liouville-von Neumann 
type (with a Hamiltonian ±iK/2) : 


dAi) 

du 


I [!<. .4131 


dvA 

du 


1 

2 


[K , vA] 


(2.34) 


Note however that the associated boundary conditions on AV and VA cannot 
be expressed in a simple form. Indeed the conditions (|2.5| ) and ( p.ll| ) do 
not involve V and A evaluated at the same argument, but V{0) and A{P) 
evaluated at the initial and hnal limits of the integration interval. 

In Sects. 3 and 4, A{u) and V{u) will be chosen within the trial class of 
exponentials of quadratic forms of creation and annihilation operators. All 
the properties of the variations 6A and dV discussed above are then satished; 
in particular, the quantities 5AA~^, A~^ 6A, 6VV~^ or 'D~^ 5V span the set 
of quadratic forms. 


2.3 Cumulants of Conserved Observables 

The variational principle (p.4|) provides in general rather complicated expres¬ 
sions for the correlations, and more generally for the higher order cumulants 
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of the observables Q^. In this Subsection we consider a single observable Q 
which is conserved, that is, which commutes with the operator K. (The ex¬ 
tension to several observables commuting with K, but not necessarily among 
themselves, is readily performed by replacing by J2-y^-yQ'y) In fhe exact 
case, one has the obvious relation 


MO = IV, e-iQ = IVi 


(2.35) 


which tells that the exponential of the characteristic function o associated 
with the observable Q is equal to the partition function corresponding to the 
shifted operator 

K\^)^K + ^Q/P . (2.36) 

We shall see now that the stationary solution of our variational principle 
satishes also the property (p.35|) when the trial spaces for the operators A 
and V are invariant under left and right multiplications by exp(—AQ), where 
A is any c-number. This invariance property is, for instance, satished when Q 
is a single-quasi-particle operator and when A and V are, as in Sects. 3 and 
4, exponentials of a quadratic form. We shall also encounter the property 
( |2.35|) in a different context in Sects. 5 and 6. 

Let us indeed consider a trial set {A{u), V^u)} satisfying the boundary 
conditions and ( |2.11| ) with A = exp{—^Q). We note that the two 


operators A'{u), V’{u) dehned by 

A'{u) = A{u) e<^/2/^ , V'{u) = V{u) e-<^/2/^ , 

(2.37) 

obey the boundary conditions ( p.5|) and (|2.11|) with A = 1. Moreover the 
value of the functional (1^.41) calculated with A'{u), V'iu), A = 1 and K' is 
equal to the value calculated with A{u), A = exp{—^Q) and K. By 

sweeping over the trial set {A{u), T>{u)} and using the above invariance of 
the variational spaces, we obtain 

exp (p{e-^^, /3, k} = Js{e-^^, /3, K} = J,{i, A, K + ^Q/A} , (2.38) 
which, within our approach, is the variational counterpart of the relation 


(|2.35|) . As a consequence, for any trial space satisfying this invariance, all 
the cumulants at temperature 1//3 of a conserved observable can be obtained 
from a calculation of the partition function for the shifted operator K' dehned 
In particular, the mean value {Q) and the variance AQ are given 


m 
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by 


(Q) = 


^(^{e f3, k} 


= ft A'+ AQ} 

g=0 


A=0 


AQ^ = 




^ip{e-^Q,P, K} 


«=o 


1 d^ 


if{i,(3,k + XQ} 


A=0 


(2.39) 

which are consistent with exact thermodynamical relations 

An example of special interest concerns the particle-number operator N 
in the grand-canonical equilibrium {K = H — yiN). To obtain the cumulants 
of N (i.e. A = exp(—^iV)) it suffices to consider the partition function in 
which fi is replaced by /i — A with \ = ^/j3. Then Eqs. ( p.39|) become 


[N) = 13, H - ^N] 




(2.40) 


These equations, valid in any variational space that satishes the aforemen¬ 
tioned invariance, coincide with the usual thermodynamic relations for the 
expectation value and the fluctuations of the particle-number operator. Sim¬ 
ilar relations arise for any conserved observable such as the momentum or 
the angular momentum. 


2.4 Connection with the Standard Variational Princi¬ 
ple 

While for 'D{u) the exact solution exp(—uA') has a trivial, exponential de¬ 
pendence on M, the exact solution (|2.12|) for A{u) depends on m in a more 
complicated fashion since lnyl,(M) is a linear function of u only if the operators 
A and K commute, as was the case in Sect. 2.3. For restricted variational 
spaces, due to the coupling between the differential equations ( p.6D and (|2.8|) , 
we expect in general both the approximate solutions T>{u) and A{u) to dis¬ 
play a non-trivial n-dependence. This is actually the feature which provides 
its flexibility to the method. 

Nevertheless, for = 0 (A = 1), that is, when we only want a varia¬ 
tional evaluation of the thermodynamic potential </3(0), a simplihcation oc¬ 
curs whenever the trial spaces satisfy the following properties (which will 
hold in all the Sections 3 to 6 and Appendix A) : (i) the trial sets for the 
operators V{u) and A{u) are identical; (ii) if V belongs to this set, so does 
cP" where c and a are any positive constants. (The power a of the operator 
t) is well-defined in a basis where it is diagonal, even when it has vanish¬ 
ing eigenvalues. This remark will allow us to use the results of the present 
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Section not only for independent-quasi-particle density operators (Sect. 4.1) 
but also in Sects. 5.4, 6.2 and 6.3 where we will consider projections of these 
density operators onto subspaces of the Fock space.) The properties above 


enable us to make the following Ansatz, suggested by Eq. (p.l2|) for A = 1 : 


Vo{u) = e 


-uHn 


Ao{u 


^ -{(3-u)Ho 


(2.41) 


where the trial operator Hq does not depend on u. As a consequence, we 
have 


Ao{u)'bQ{u) = Po(m)A(m) = VM = e 


zA 


(2.42) 


where we dehned Z as the normalization Tri e The variational expres¬ 
sion (p.4D then reduces to 


X = Z(1 - InZ) - ZTri A (In A + (3k) 

By maximizing ( p.43|) with respect to Z , we obtain 

InX = InZ = -Ti'i AhiA - (3 Tri A A 


(2.43) 


(2.44) 


which shows that, under the conditions of the present section, the entropy S 
dehned by ( |2.28|) is equal to 


S = -Tri AlnA 


(2.45) 


and is related to T>{j3) by Eq.(|L2[). The standard variational principle for 
thermodynamic potentials is thereby recovered ; indeed we hnd (p(0) as the 
maximum with respect to A (under the constraint Tri A = 1) of the right- 
hand side of ( p.44|) , in agreement with the fact that InZ is the Legendre 
transform of the entropy. This reduction, however, takes place only for A = 1 
or = 0. 

The analysis of the previous Subsection showed that for an operator Q 
commuting with K, the characteristic function (p(^) is equal to the thermo¬ 
dynamic potential associated with the shifted operator K' = K+^Q/(3 when 
the trial spaces satisfy the invariance property stated in Sect. 2.3. In this 
case the cumulants of any order, and in particular the correlations, can be 
deduced from the standard variational principle (|2.44|) with K replaced by 

k'. 
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2.5 Variational Principle for ^(w) Alone 


The Ansatz ( |2.41|) , which restricts the parametrization of the two M-dependent 
variational quantities T>{u) and A{u) to a single operator Hq, is not suited 
to the evaluation of the characteristic function (p(^) for 7 ^ 0 (and therefore 
of cumulants of non-conserved quantities). In the following Sections we shall 
therefore exploit the freedom allowed by the full set of the variables 'D{u) and 
A{u). Nevertheless the form of the exact solution ( | 2 . 12 |) for A{u) suggests 
another possibility : instead of letting A{u) and T’(m) vary independently, 
one can constrain them according to 


A{u) = - u)AV^/^{(3 - u) 


(2.46) 


Inserting this trial form into the functional ( |2.4|) leads to the variational 
expression 


Tn AV{(3) - r dw Tn - u)AV^/\p - u) 

Jo \ dti 


(2.47) 


+ - [kV{u) +V{u)k] 


in terms of a single variational quantity T>(u) obeying the boundary condition 
(|2.5| ). Had we not symmetrized the Bloch equation, we would have found 
the alternative functional 


Tn Av{p) - / du Tn Av{p - u) 


'dVju) 

du 


+ kV{u) 


(2.48) 


whose stationary value with respect to Viu) yields the wanted quantity 
Tr Ad. a special case of this variational principle has already been pro¬ 
posed and applied in connection with the Monte Carlo method |^; it is 
obtained from ( p.48|) in the case of canonical equilibrium k = H hj taking 
for A the dyadic |i?)(i?'|, where i? is a point in the 3A^-dimensional space. 
The stationary value of (p.48|) is then the density matrix itself in the R- 
representation. 

If we further specialize T’(m) to depend exponentially on u, as does the 
exact solution, in the form 


V{u)=e--^^/A , 

both functionals (p.47| ) and ( |2.48| ) become 


Tri Ae~-^ ( i 




(2.49) 


(2.50) 
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We thus recover the variational expression (3.20) of Ref.p3 
trial quantity is the operator A4 


In (|2.50| ) the 


For unrestricted variations of A4, the sta¬ 
tionary condition is Af = (3K, and the stationary value is the characteristic 
function Tri Ae~^^. It has been shown in that, unlike the previous vari¬ 
ational functionals, ( 2.50 ) is maximum for Ji4 = [3K under rather general 
conditions. 
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3 Extended Finite-Temperature HFB Approx¬ 
imation for Characteristic Functions 

From now on we will consider systems of interacting fermions for which pair¬ 
ing effects are important. In Sect. 3.2 we therefore apply the variational 
principle of Sect. 2 to derive consistent coupled equations which are adapted 
to the evaluation of the characteristic function ( 0 ) and which account for 
pairing. We deal, in this Section and Sect. 4, with grand canonical equi¬ 
librium, letting K = H — piV and taking traces in the Fock space. The 
formalism will also apply if K includes additional constraints such as —cuJ, 
where J is the ^-component of the angular momentum. More generally, for 
a hnite system, one may want to control its position in space by including in 
the set Clk the center-of-mass operator, or, if it is deformed, its orientation 
by the quadrupole-moment tensor. We first introduce some notation and 
dehnitions which will be used throughout the present article. 


3.1 Generalities and Notations 

For the sake of conciseness, it is convenient to introduce 2?7,-dimensional 
column and line vectors dehned as 


/ cq \ 


a = 




\ ) 


= (^a\ ■ ■ ■ a^n ai ■ ■ ■ a. 


.) . (3.1) 


where a^i and a* are the creation and annihilation operators associated with 
a hxed single-particle basis (i = 1, n). If we denote by ax the A-th component 
{A = 1, 2n} of the column vector a, the fermionic anticommutation rules are 
given by either one of the relations 


[cXA, «!']+ “ ^AA' , [«A, Q(A']+ — CTAA' 


(3.2) 


where a is the 2n x 2n matrix 


cr = 


0 1 
1 0 


(3.3) 


In what follows, as in the HFB theory, operators having the form 

T = exp(—/ — \cx} Ccx) , (3.4) 
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where I is a c-number and C is & 2n x 2n matrix, play a crucial role. As a 
consequence of the anticommutation relation (|3.2| ), we can restrict ourselves 
to matrices C which satisfy the relation 

aCa = , (3.5) 

where denotes the transposed matrix of C. Instead of I and C one can 
alternatively use the c-number r and the 2n x 2n matrix T defined by 

T = e~^ , T = e~^ , (3.6) 


to specify the operator (|3.4|) . Another equivalent parametrization involves 
the normalization (, 

C = Tr T = e~^ ^det(le“'^)^ = r exp ||tr 2 „ ln(l-h T) | , (3.7) 

and the generalized reduced density matrix 77 (or matrix of the contractions), 
Tr axTay f 1 \ 


77aa' = 


Tr T 


e'^ -I- 1, 


AA' 


T 

r + 1 


XX' 


(3.8) 


In Eq. (^), the symbol tT 2 n denotes a trace in the space of 2n x 2n matrices 
while Tr is the trace in the Fock space. The property (|3.5|) implies the 
relations : 


aT-V = T' 


aTZa = 1 — 71' 


(3.9) 


As a consequence, the 2n x 2n matrix 77 can be written in the form 


77 


P \ 

K+ 1- p'' ) 


(3.10) 


where the nxn matrices p, K- and /?+ correspond to the normal and abnormal 
contractions : 

f pij = Tr aifa^j/Ti f , 

< = Tr ttifaj/Tr f , (3.11) 

[ t^+ij = Tr a^jTa't'j/Tr T . 


The matrices K- and k+ are antisymmetric. When the operator T is hermi- 
tian, one has p = p\ and the matrix 71 = TV reduces to the usual 

HFB form. 

Operators of the type (p.4|) allow the use of a generalized Wick theo- 

They also form a non-abelian 


rem, 


even when they are not hermitian 


multiplicative group, the algebra of which is characterized by the relations 
(|3.2(]| - |T^) . These simple algebraic properties (more details can be found in 
Z7| and the Appendix of [^) will be used extensively below. 
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3.2 The Variational Ansatze and the Associated Func¬ 
tional 

Let us come back to our problem. We assume that the observables are 
operators of the single-quasi-particle type (this restriction is not essential and 
can be removed at the cost of some formal complications (see pi||)): 


Q'y = q-y + \CX^Q^CX 


(3.12) 


so that the operator A(^) (see Eq. ( |1.6|) ) belongs to the class (|3.4|) . We thus 
write it as 

A(,^) = = exp(—, (3.13) 


with 




(3.14) 


Now we choose the trial spaces for the operators 'I^('u) and A{u) entering 
the functional (|2.4|) . This choice determines our variational approximation. 
We take as trial sets : 


Vin) = fM(u) , A(u) = f l“l(u) 


(3.15) 


with 


d’^\u) = exp(—/^'^'( m) — d'^{u)oL) , (3.16) 

f [“!(«) = exp(-/[“l(M) - \(dd°^{u)cx) . (3.17) 

The 2n x 2n matrices d^{u) and (constrained to satisfy the relation 

(|3.5| )) and the c-numbers /[‘^^(m) and /[“^(m) constitute a set of independent 
variational parameters. Using Eqs. (HO. we define the coefficients 

the normalization factors ^[“1, the matrices T[“l, and the reduced 
density matrices 77[“1, associated with the operators T[“l and 
The boundary condition (|2.5|) on 77(0) implies the relations: 

(3.18) 


ll'4(0) = 0, £'''’(0)=0 


or equivalently 

r['^'(0) = l, r['^l(0) = l 


(['^1(0) = 2”, 77['^1(0) = - . (3.19) 


The products P(m).4.(m) and Ai^dV^u) occur in the functionals ( p.4|) and 
(|2.7| ). Taking advantage of the group properties of the operators ( |3.4|) , 
we write 

AV = , VA = f [^“1 , 


(3.20) 






with 


f [“'^1 = exp(-/[“'^l - ^cx^&’^cx) , (3.21) 

f , (3.22) 

where the M-dependence of the c-numbers and and of the matrices 
£[ad] £[da] pggj^ indicated; nsing the notation ( p.6|) , one has the 


relations 




= = e' 


./H _ i[d\ ^ 


_£[ad\ 


(3.23) 

(3.24) 

According to Eqs. (|3.7|) and (|3.23| - ^)24D , the normalization of VA and AT), 
which we denote in shorthand by Z, is given at each “time” u by 


_ j-[ad\ _ j-[a\j-[(]\ 


_£[da\ _ rj'^da] _ rj'[d\t-j-[a 


Z =1)1 AV = =1)1 VA = 

= e-('‘“'+'''^)exp|itr,„ln(l + e-£'“ 

exp{|tr 2 n ln(l + rHrl'^1) | 


.£^1. 


(3.25) 


= 


[d] 


As we already know from Sect. 2.2, the normalization Z does not depend 
on u at the stationary point and it is a variational approximation for the 
qnantity exp(p that we are looking for (see below Eqs. (|3.48| - |3.50|) ). It thns 
plays the role of a generalized partition fnnction which rednces to the grand 
canonical partition fnnction when the sonrces are set to zero. 

We can now write the fnnctional (p.4|) (which provides the characteristic 
fnnction ( p.l|) ) in terms of the variational parameters T[“l. The 

term TiAdV/du is obtained by differentiating and in Eq. ( 3.25 ) with 
respect to u, while holding A°^ and hxed, which yields 


Tr A^ = Z 
au 


—- + -tT2n ' Ti^ ' - 

du 2 du 


(3.26) 


The last two terms in the integrand of the fnnctional ( |2.4|) are snpplied by 
a straightforward application of the generalized Wick theorem. They involve 
the operator K = H—fiN where the Hamiltonian H of the system is assnmed 
to have the form 


and where 


H ^ ) ^ij Aidj T ^ ^ ) ^ijkl AiAjOjlOj}^ 

ij=l ijkl=l 


N = Y^ Aidi 

i=l 


(3.27) 


(3.28) 


29 













The matrix elements Bij of the one-body part of H and the antisymmetrized 
matrix elements Vijki of the two-body interaction satisfy the relations 

Bij = , Vijki = —Vjiki = —Vijik = . (3.29) 

As mentioned at the beginning of this Section, the one-body part associated 
with the matrix B may include, in addition to the kinetic energy, some other 
single-particle operator whose average value is kept hxed. One hnds 

Tr KAV = , Tr KVA = , (3.30) 

with 

E{JV\ — ^(0 {Bij — n^ij) Pji ^(0 Vijki {‘^pkiPij ~ ^+ij^^-ki) ■ (3.31) 

ij ijkl 


When p = p^ and = kA, the expression ( |3.31|) is formally identical to 
the HFB energy. However, the quantities E{TZ^°‘^} and are not 

necessarily real since and do not have to be hermitian. 

Altogether, for the trial choices (|3.16|) and (|3.17|) , the functional 
reads 


I{A, V} = 


f dlnrl'^1 

du Z 


10 \ dti 

1 -1 dT^'^l 1 ' 

2 dri 2 


where use has been made of the relation 

According to Eq. ( 3.25|) , the boundary term is given by 

ZA = Tr A{i)VA) = Tr fV\VA) 

= jAAid] exp {itr 2 „ln(l + r[^lr['^](;5)) } 


(3.32) 

(3.33) 

(3.34) 


where the c-number and the matrix which specify the characteristic 
function we are looking for, are expressed in terms of the sources A tiy 
Eqs. (|3.12H3.1^) and ( ^^ . Similarly the form (^^ of J reads 


X{A, V} = - ZA) + Z{d) + du 

Jo 


/dlnrl" 

du 


1 dTb] 1 

--tT2n ^[“1 —- 

2 dri 2 


(3.35) 
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In the functionals (|3.32|) and (|3.35|) , the n-dependence has been omitted in the 
variational parameters and r [a]^ ^[a] 

, as well as in the normalization 
Z given by Eq. (|3.25D and in the contractions associated with 

respectively; the latter quantities are related to and 
through (|3.24|) . 


3.3 The Coupled Equations 

The next step is to write the equations expressing the stationarity of the 
functional ( fl.32| ) or ( |3.35 ). We derive them directly from these expressions, 
but could also have used the general equations (|2.6|) and (^^p.9|) . 

Requiring the stationarity of ( p.32|) with respect to (in the func¬ 

tional (|3.32|) , r[“](«) only appears as a proportionality factor in Z{u)) gives 
us the equation of evolution for = —lnr[‘^(M) : 


^ = -tr2„r[‘'l 

dn 2 dn 2 


(3.36) 


which we shall solve explicitly below (Eq. (|3.46| )). 

We now require the stationarity of the functional (|3.32|) with respect to 
T[“](m), which appears through and Z{u). Eq. (|3.36|) implies 

that we can disregard in this calculation the variation of Z{u). If we note 
that the variations of (m) and (u) resulting from those of T[“1 (m) 
satisfy ^ as a consequence of ( p.33| ), we readily hnd 

the differential equation 


where T-[{TZ} has the same formal definition, 

6E{n} = ltT2nn{n}6n , 


(3.37) 


(3.38) 


as the HFB Hamiltonian, although 7l{7l} is not necessarily hermitian. For 
the Hamiltonian (|3.27|), the 2n x 2n matrix 71 takes the form 


= ( t -h- 


(3.39) 


with 


hij — 
= 

^+ij - 


{^Bij p,(fjj) 7 '22ikl ^ikjlPlk 
2 Efc/ ^ijklf^-kl 1 
2 Thkl ^+kl^klij 


(3.40) 
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As a consequence of Eq. ( p.29| ), the matrices A_ et A_|_ are antisymmetric. 

(3.41) 


This entails for TL the same relation as ( p.5|) : 

aTia = —07 . 


which is consistent with a57la = —57V- and Eq. (|3.38|) . 


Likewise, the stationarity conditions of the functional (|3.35|) with respect 
to and V^{u) provide, for 0 < u < /3, the differential equations 


^ = -tr2nr[“l 
du / du 

for = —lnr[“l(-u) and 




ut'H 1 

= - [Lf}r[“1 + }] , 

0 . 7 / 2 


(3.43) 


for T[“1 (m). 

Finally, by rendering the functional (|3.35|) stationary with respect to 
r['^l(/3) and r['^l(/3) (see Eqs. ( p.l3|) and (|3.6| )), one obtains the boundary 
conditions for A{u) : 




(3.44) 


The main equations to be solved are the differential equations (|3.37| ) 
and (|3.43|) for V'^{u) and T[“l {u) which are coupled through the matrices 
-^{-^[ad]} 7i{7V°^}. Moreover one has to deal with the mixed boundary 

conditions 

T['^1(0) = 1 , r[“l(/3) = . (3.45) 

All our variational quantities depend upon the sources (see Eqs. (^.12[- 


3.14|) and (|3.6| )) through the boundary conditions (|3.44|) . 

One notes that the evolution equation ( ^.43]) of can be deduced from 
Eq. ( p.37| ) for by an overall change in sign of the right hand side and 
by the exchange a ^ d, as could have been inferred from the comparison 
between the forms (|3.32|) and (p.35|) of the functional I. (The same is true 


for the equations (|3.36|) and (|3.42|) governing and if one takes account 
of Eqs. ( p.37|) and ( 3.43| ).) 

Once Eqs. (p.37|) and (|3.43|) are solved, the c-numbers and V^{u) 

are obtained by integration of Eqs. (p.36| ) and ( |3.42|) whose r.h.s. depend 
only on and . Indeed, using Eqs. (p.37|) and ( p.43|) with the boundary 
conditions /['^(O) = 0 and we find 


lW{u) = p du {-|tr 2 n 


(3.46) 
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/[“I (m) = {(3) - (m) . (3.47) 

From the equations for T[“l, /[“I, together with the dehnition 

(|3.25|) and the relation (|3.38|) , one can readily verify the conservation laws 


(|2.17|) and (|2.15|) which now read, respectively : 


dn 


= 0 


dn 




(3.48) 


One can also verify that the integrands of the functionals (|3.32 ) or ( 3.35 ) 
vanish. As was shown in Sect. 2, these properties directly result from the 
feature that variations 5A oc A and 5T> oc T) are allowed in the trial classes 
dSA^^). 

We are looking eventually for the characteristic function (p{A(^),/?, 77} 
which is given by the logarithm of the stationary value of the functional 
(|3.32|). This value is provided by the boundary term dehned in ( B.34| ). 


In this formula, the matrix and the c-number rt‘^(/3) = exp [—/I 

are furnished by the solution of the coupled Eqs. (|3.37|) , (|3.43|) and by (p.46|) , 
respectively. Since, according to Eq. (|3.48|) , the normalization Z (dehned in 
(|3.25|) ) does not depend on u we have the relation 


<7(0 = = InZt^l = \nZ{u) , 0 < u < (3 


(3.49) 


In particular, Z can be evaluated at n = 0 (see Eq. (|2.19|) ); from Eqs. (|3.25|) 
and (|3.19|) we hnd 


InZ = -/[“l(o) + kr 2 „ln(l + r[“l(0)) 


(3.50) 


where /H(o) = /Al + /F](/3) is explicitly given by Eqs. (|3:T5H^ ) once the 
coupled equations (|3.37|) and (|3.43|) which yield T[“](0) have been solved. 

Let us add a few remarks about the properties of these last equations. 
In the derivations of this Section (as well as in Sect. 2), the hermiticity of 
the operators V^u) and A{u), as dehned in Eqs. ( |3.15| - |3.l7| ), was nowhere 
assumed. If, for some value of u, the operators P and A are hermitian, 
i.e. if and = A°^*, it follows that 

and hence that . When the operator A 

is hermitian, the inspection of Eqs. (|3.37|) , (|3.43|) and ( |3.45| - p.47|) shows that 
these equations preserve the hermiticity of and and the reality of 
A^^ and A^\ In this case, the stationary value Z given by (|3.50|) is also real. 
If A is not hermitian, as would happen for the characteristic functions with 
imaginary ^’s, the boundary condition ( p.44|) prohibits A{u) and T>{u) from 
being hermitian. 

Combining Eqs. (p.37|) and (|3.43|) , one obtains the evolution equations 
for the products and , and hence for the 
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HFB-like contraction matrices and associated (through Eq. ( |3.8| )) 
with the operators and The calculation relies on the property 

57^ = (1 - 7^)(5r(l - 7^) , (3.51) 


and it leads to 




= . (3.52) 


These equations, the single-quasi-particle reductions of Eqs. ( p.34|) , involve 
commutators. They are reminiscent of the time-dependent-Hartree-Fock- 
Bogoliubov equations (see Sect. 5 of the Appendix A), iiw replacing the 
time t. However, here and are in general not hermitian 

and u is real. A consequence of the simple form (p.52|) is that the eigenvalues 
of and remain constant with u, in contrast to those of and 
It also follows from Eqs. ( p.52| ) and (p.38|) that each of the quantities 
and is conserved, a hner result than Eq. (|3.48|) . In spite of 

their form, one should not infer from Eqs. (p.52|) that the matrices 7^1“'^! and 

^[da] 

are decoupled. Actually the constraints (p.l9|) and (|3.44|) do not entail 
any explicit boundary conditions for 7^1“'^! and 7^1'^“!. Moreover, Eqs. ( 3.52 ) 
are not equivalent to Eqs. (p.37|) and (p.43|) since, except in very particular 


cases, 7^[“] and are not determined uniquely by knowledge of and 

7^ [da] _ 
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4 Expansion of the Extended HFB Approxi¬ 
mation : Fluctuations and Correlations 


In this Section we derive the thermodynamic quantities, the expectation 
value of the observables, their fluctuations and the correlations between these 
observables from the characteristic function as given by either one of the 
formulae (|3.34|) or (|3.50|) . We shall see that this evaluation does not in 
fact require the full solution of the coupled equations (|3.37|), (|3.43|), 


and ( p.47| ) for hnite values of the sources Indeed these quantities can be 
obtained through a prior expansion of the coupled equations in powers of the 
sources followed by a solution of the resulting expanded equations order by 
order. This procedure turns out to overcome the difficulties associated with 
the mixed boundary conditions. As expected, the usual HFB equations are 
recovered in the case A = 1 , when the sources vanish (Sect. 4.1.1). These 
equations are also relevant in the calculation of the average values of single- 
quasi-particle operators (Sect. 4.1.2). In Sects 4.2 and 4.3 we work out the 
second order in the sources ; this yields an explicit variational expression for 
the fluctuations and the correlations which differs from the naive application 
of the HFB approximation. In the present Section, the increasing orders of 
the expansion with respect to the will be indicated by lower indices (0, 
1, ...) assigned to the considered quantities. 


4.1 The HFB Approximation Recovered 

4.1.1 Thermodynamic Quantities 

To zeroth order in the sources (i.e., A(^) = 1) our problem amounts to 
evaluating the partition function Z = Tr I?o(/3). We expect that the coupled 
equations of Sect. 3.3 will yield the standard HFB approximation for Z, 
since the trial spaces for V^u) and A{u) satisfy the conditions required in 
the discussion of Sect. 2.4. Let us verify this point. If we associate with the 
contraction matrix TZq the single-particle Hamiltonian Tio dehned by 


the usual optimum HFB matrix TZq is given, at temperature 1//3, by the 
self-consistent equation 

Ho = Lf{7^o} , (4.2) 


where H{TZ} was dehned in Eqs. (|3.39| - |3.40|) . Because the relation (^) 
entails the commutation of 7-f{7^o} with TZq, Eqs. (|3.52|) are satished by the 
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M-independent solution 


= n^t\u) = 7^o 


(4.3) 


In the equations (|3.37|) and ( p.43|) for Tq^{u) and Tq°^{u), the effective Hamil¬ 
tonians do not depend on u since } = Tio- Moreover, the 

boundary conditions ( |3.45|) reduce to = 1, Tq^\(3) = = 1. We 

thus hnd : 


TP{u) = e““^o , . (4.4) 


Eqs. (|4.3H4.4|) are the single-quasi-particle reductions of the general relations 
(|2.42|) and (p.41|) . One also notes that TZo is equal to 7l^^{/3), the contraction 
matrix associated with the state Vo{(3). 

Using the solution ([4.4|) of the coupled equations, we can now calculate 
/[“](0) from Eqs. and hence InZ from Eq. p.SOj ). The integrations 

over M, which can be performed explicitly, yield InZ from which according to 
Eqs. (|2.27|) and (p.28|) we dehne the thermodynamic grand potential Fq and 
the entropy S'{7^o} 


InZ = S{7Zo} - /3E{7^o} = -/3Fg{7Zo} 
The entropy is given by 


(4.5) 


S{7J} = -ltr. 2 „ [TllnR + (1 - K)lii(l - K)] . (4.6) 

It is thus the entropy of a quasi-particle gas characterized by the contraction 
matrix TZ. Owing to the symmetry relation (|3.9|) , the two terms under the 
trace in (|4.6|) are equal. The HFB energy E{TIq} = {H — /iiV) is given by 
the expression (t3.31|) . 

As expected from the discussion of Sect. 2.4, Eq. ( |4.5| ) coincides with the 
outcome of the variational principle (|2.44|) , which reads here 

InZ = InTr Vo{f3) = [S{n} - /3E{7^}] = -f3 Min FciTZ} . (4.7) 


The stationarity condition of (|4.7|) is obtained from (|3.38|) and from the hrst- 
order variation 

6S{n} = ^tT 2 n In ^ ^ 5n (4.8) 

of the reduced entropy (|4.6|) ; it thus gives back the usual self-consistency 
condition (T^). Therefore, at zeroth order in the sources, we have recovered 
the standard self-consistent HFB approximation for the thermodynamic po¬ 
tential. 
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4.1.2 Expectation Values of Observables 


The variational approximation for the expectation values (Q 7 ) in grand 
canonical equilibrium is in principle given by the hrst-order terms in the 
sources ^ 7 . However, the knowledge of the zeroth-order solution H). (B 
is here sufficient since, according to Eq. ( p.23| ) and to the definitions (|3.12|) 
and (|T^), we have 


{Q')) ^ ^ 7^0 

Tr Vq[P) 2 


(4.9) 


For the expectation value of the operator K, we can apply Eq. (|2.26|) 
which was obtained by derivation with respect to the parameter /?. Using 
Eqs. (13.301) at zeroth order and ( |4.3|) , we find 


(k) = £{R„} = -UlnZ 


(4,10) 


where, as in Eq. ([4.9|) , the expectation value is taken over Vq^P). 

For the calculation of (N), we can use either ( |4.9|) or the discussion in 
Sects. 2.2 and 2.3. From Eq. (p.31|) (or (|2.4(]|) ) we get 


« fi \ \ d 

(N) = tr„p = - + 


(4,11) 


where, as in Eq. (p.l2|) , the operator N can be represented by the c-number 
n /2 and the 2 n x 2n matrix 


M = 


1 

0 


(4,12) 


The formulas ( |4.9| - |4-11]) are those given by the standard approximation, 
which consists in evaluating expectation values in the HFB state Po(/^)- Our 
variational procedure thus does not yield anything new when the character¬ 
istic function is expanded to first order in the sources. This is not the case 
at the next order. 


4.2 Fluctuations and Correlations of Conserved Ob¬ 
servables 

4.2.1 The Particle Number Operator 

As a first example of a conserved observable (i.e., which commutes with K), 

^ ^ 2 

we consider the particle number N and we evaluate its variance AiV . A 
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naive method would consist in taking the difference between the expectation 
value {N"^) and the square both evaluated by means of the Wick’s 

theorem applied to the state Vq{P). This would yield 

(7V2) - (iV)2 = - 7^o)A^7^o, (4.13) 

where TZo is the contraction matrix (|4.1| - |4.2| ) corresponding to Vq{(3). 

However, the state T>q{13) is variationally fitted to the evaluation of the 

^ 2 

thermodynamic potential Fq, not to the evaluation of AN . In agreement 
with our general philosophy, we should rather rely on a variational princi¬ 
ple suited to the evaluation of the characteristic function ip. The variational 

- 2 

approximation for the fluctuation AN is then supplied by the second deriva¬ 
tive of (p with respect to the source associated with the observable N. The 
result is expected to differ from ([4.13| ), since the optimal values of the trial 
operators T’(m) and A{u) are thus dependent on the sources 

Actually, under conditions on the trial classes that are obviously satisfied 
here, we have already determined in Sect. 2.3 the result of this evaluation : 


-2 ld{N) Id^Fg 

P dp P dp^ 


(4.14) 


where Fq and {N) are the outcomes of the variational calculation at ze¬ 
roth and first orders in the sources, given by (|4.5|) and (14.11|), respectively. 

^\ 2 

The identities ( [4.14| ), satisfied by our approximation for AN , agree with 
exact thermodynamic relations which are violated by the naive result ([4.13|) . 
Inserting ([4.11j) into (^4.14|), we have 


1 ^.dTZo 


(4.15) 


The derivative dTZo/dp is supplied by Eq. (^]2|) which defines TZq self-consistently; 
the parameter p enters this equation through the contribution —pM to 
7-f{7?.o}. The expression (|4.15|) reduces to (|4.13|) only for a system of non¬ 
interacting particles, with [A/", T^o] = 0, and is otherwise non trivial due to 
self-consistency. We shall write it explicitly below (Eq. (|4.27|) ) in a more 
general context. 

The result ([4.15|) arising from of our variational treatment is more satis¬ 


factory than the naive expression (|4.13|) . We noted that it meets the thermo¬ 
dynamic identities ([4.14|). Moreover, for a system whose pairing correlations 


do not disappear at zero temperature, Eq. ([4.13|) has a spurious finite limit as 
P goes to oo, whereas the variational estimate (|4.15|) for the particle-number 
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dispersion vanishes in this limit as 1 //?, with a coefficient eqnal to the density 
of single-particle states at the Fermi level. It is therefore consistent with the 
expectation that, for hnite systems with or withont pairing, resnlts from the 
grand-canonical and canonical ensembles mnst coincide at zero temperatnre. 


4.2.2 Characteristic Function for Conserved Single-Quasi-Particle 
Observables 

The ideas above are readily extended to correlations between any pair of 
single-qnasi-particle observables and Qs commnting with the operator 
K. Indeed, it follows from the discnssion of Sect. 2.3 that the evaluation of 
the characteristic function, and therefore of the cumulants, can be reduced in 
this case to the calculation of a thermodynamic potential if the trial spaces 


for A and T> are invariant with respect to the transformation ( p.37|) . Within 
the spaces ( p.l 6 |) and ( p.l7| ) considered in this Section, this invariance is 
satished for operators of the type (|3.12|) . 

The identity ( |2.38|) then shows us that the variational evaluation of the 
characteristic function (p(^) amounts to the variational evaluation of the par¬ 
tition function associated with the shifted Hamiltonian K' = 'C 7 Q 7 / 

The calculation is thus the same as in Sect. 4.1.1, except for the change of 


K into K'^ and Eq. (|4.7|) replaced by 


V9(0 = Max7^[5{7^} - PE{n} -Y.U<h + Q,n) ] 


(4.16) 


The optimum matrix IZ' is characterized by the extremum condition (|4.16|) , 
which yields the self-consistent HFB equations 


= 




+ 1 




(4.17) 


The characteristic function ([4.16|) depends on the ^’s both directly and 


through the matrix TZ'. From its hrst derivatives we recover the expectation 
values (^]9|). The second derivatives, taken at .^ = 0, yield the correlations 
between the Q.y’s. In order to write their explicit expression, setting IZ' = 
TZq+51Z, we shall expand S{7Z'}, E{TZ'} and Eq{TZ'} around TZq up to second 
order. To this aim, we introduce a condensed notation. 

Up to now, we have regarded TZxy, Tiw as well as other quantities de¬ 
noted by a script capital (Q, £, T, 57Z, Hq, etc.) as 2n x 2n matrices. We 
hnd it convenient to consider them alternatively as vectors in the Liouville 
space, with the pair (AA') playing the role of a single index (see Appendix, 
Sect. A.l). Thus, the hrst-order variations ( |4.8| ) of S{TZ} and ( |3.38| ) of E{7Z} 
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appear now as scalar products that we shall denote by the colon sign : which 
stands for |tr 2 n . Likewise the second-order variations will generate matrices 
in the Liouville space, with two pairs of indices (A A'). 

With this notation we can write, for IZ' = TZq + 671, the expansions of 
the HFB entropy S'{7^}, energy E{71} and grand potential Fq{71}, dehned 
by Eqs. ( |4.6|) , (|3.31|) and (|4.5|) , respectively, as (see Appendix, Sect. A. 2 ) : 


~ 5{7^o} + f3no : 6n + \5TZ : So : 57^ , 
where /3?-fo stands for ln(l — 1Zq)/71q as in ( [4.1| ) ; 

E{7i'} ~ E{7^o} + L^{7^o} : 511 + \5n ; Eq : 57^ , 
where l-i{7lo] is defined by ( p.38| ) ; 

FG{7^'} ~ FG{7^o} + \5n ; Fo : 57^ , 
where we made use of the stationary condition (|4.2|) and where 


Fq = Eq — —So 


(4.18) 

(4.19) 

(4.20) 

(4.21) 


In Eqs. ( [4.18| - |4.20|) , the second derivatives So, Eo and Fo, as well as the first 
derivatives [37io and 7-f{7?.o}, are functions of the HFB density TZq. In the 
A-basis of Sect. 3.1, So, Eo and Fo are 4-index tensors ; we regard them here 
as (symmetric) matrices in the Liouville space (AA'). The explicit expression 
of So(AA')(A(/i') is given in Eq. ( |A-1(]|) . Contraction of such a matrix with a 
vector proceeds as 


(Sq : 671)xy = ^ SofAApfw') • (4.22) 

fj.fl' 


In the Appendix A (Sect. A.4), we endow —Sq with a geometric interpre¬ 
tation. It appears naturally as a Riemannian metric tensor which defines a 
distance between two neighbouring states characterized by the contraction 
matrices TZq and TIq + 671. 

The correlation C.ys between the two observables Q7 and Qs (defined as 
in Eqs. (|1.7| - |1.8|) ) is obtained from the expansion of the expression (|4.16|) 
around TIq : 

(p(0 ^ InZ -YjUQi) - MmsnilPSn : Eq : ^ ^767 ^ • (4-23) 

7 7 

The extremum condition, 

PFo:6'JZ = -Yj^'yQ'y ^ ( 4 - 24 ) 

7 
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determines 6TZ as function of the (^’s. From Eq. (^4.23|), one has 


C^s = 


d^if 


d^s 


= -Q. 


?=o 


d^s 


5n 


(4.25) 


5=0 


Inserting in Eq. ( [4.25| ) the formal solution of Eq. ( 4.24 ), 


d 

9^5 


= ; Q, 


5=0 




one obtains 


C'yS — — Q 


/9 


Fo' 


: 


(4.26) 


(4.27) 


We shall discuss the problem of the vanishing eigenvalues of Fq in Sect. 4.2.3. 
The formula (|4.27|) encompasses the fluctuations AQ^, obtained for 7 = 5. 
Remark : The above results have been obtained indirectly by using the iden¬ 
tity (p.38|). We could also have found them directly by solving the coupled 


equations of Sect. 3.3. This is feasible, but not straightforward. Indeed, hav¬ 
ing defined IZ' and by Eq. ([4.17|) , we can first check that the solution 

of Eqs. (p.52|) is given by 

77^1 (u) = /213 /213 ^ /2p /2j3 _ 

(4.28) 

The proof uses the identity (^4.32| ) below where U is replaced by expu/lt"^!. 
We can then solve the coupled equations (|3.37|) and (|3.43|) , together with the 
required boundary conditions T['^1(0) = 1 and T^°^{(3) = exp(—which 
results in 


1213 12(3 

^^-uC^^^/2(3^-{(3-u)n'{R!}^-uh^^/2(3 


(4.29) 


We hnally recover the characteristic function (|4.16|) by means of (|3.46H3T47|) 
and ( |3.49| - P75(]|) . If commutes with TZ\ and therefore with the 

expressions ( |4.28D and ( |4.29|) simplify. Otherwise, when some (^.^-invariances 
are broken, the u-dependence in the solution ( 4.29 ) of the evolution equations 
(|3.37|) and (|3.43|) is no longer as simple as it is for ^ = 0 or for the exact 


solution [given by T>(u) = e and by Eq. (p.l2|) for A{u)]. In particular 


u 


lnT['^l(M) is not proportional here to u, but contains contributions in u^, 

. Although, in the case of commuting observables, we were able to indirectly 
solve the coupled equations (|3.37|) and (|3.43|) , it would not have been easy 
to guess a priori the complicated u-dependence of the solutions ([4.28[-^1^9D. 
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4.2.3 Broken Invariances 


Since equilibrium is described by the minimum of the grand potential Fq{TZ}, 
given by (b:-20|) , the eigenvalues of the matrix Fq dehned in Eq. (|4.21|) are non¬ 
negative. We have, however, to consider the case of vanishing eigenvalues for 
which the inversion of Eq. ([4.24|) raises a problem. A characteristic situation 
where this happens is that of broken invariances. For instance, in the mean- 
held description of pairing correlations, the particle-number N symmetry is 
broken. For convenience, in this Subsection, we shall use this operator as 
the representative of all broken symmetries, although the whole discussion 
applies to any conserved single-quasi-particle observable. 

By dehnition of a conserved observable, N satishes the commutation re¬ 
lation [K, N] = 0. When the W-invariance is not broken, N commutes with 
the approximate density operator and we have both [T^o 5 A/] = 0 and 

\Hq , M] = 0, where the matrix M has been dehned in ( 4.12|) . However, when 
the W-invariance is broken, M does not commute with TZq and Tio despite 
the commutation of N with K in the Fock space. In the HFB theory, this 
non-commutation manifests itself by the occurence of non-zero oh-diagonal 
elements in (|3.10|) and ( p.39| ). 

Let us show how the W-invariance, whether or not it is broken, is rehected 
on our variational objects. We can associate with N a continuous set of 
unitary transformations U = exp ieN which leave K invariant. They are 
represented in the 2n-dimensional A-space of Sect. 3.1 by unitary matrices 




(4.30) 


From the dehnition E{71} = Tr VK/Ti V and the commutation [K, N] = 
0 , we hnd 


E{n} = 


Tr UVU-^K 


= E{unu-^} 


(4.31) 


Tr UVU-^ 

for any TZ satisfying the symmetry condition (|3.9|) . We also have S'{7^} = 
S{U1ZU~^} since *S'{7^}, dehned by ( ^(^ , depends only on the eigenvalues 
of Ti. Thus, if Eq{JI} reaches its minimum for some TZq satisfying (|4.2|) and 
if the W-invariance is broken, the transformation hnZQU~^ generates a set of 
distinct solutions which give the same value Fq{U1ZqU~^} = Fg{7^o}- This 
implies that the second derivative Fq of Fq{TZ} has a vanishing eigenvalue. 
To hnd the associated eigenvector, we note that the hrst-order variation of 
(|4.31|) with respect to 71, together with the dehnition (|3.38|) of 7i.{7Z}, yields 


UH{7Z}u-^ = rtiuizu-^} 


(4.32) 
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Hence, letting TZ = TZo, taking an infinitesimal transformation U and nsing 


5TL = Eq : 671 (a conseqnence of the definition (|4.19|) of Eq), we find 

[Af,Ho] = Eo:[Af,no] . (4.33) 

The same argnment holds when E{71} is replaced by Fq{71}, bnt then the 
hrst derivatives of Fq vanish at 7^ = TIq and Eq. (|4.33|) is replaced by 

Fo:[A^,7^o]=0, [A^,7^o]:Fo = 0 . (4.34) 

These equations express that the commutator [Af , T^o], which does not vanish 
when the iV-invariance is broken by TZo, is a right or left eigenvector of the 
matrix Fq for the eigenvalue 0 . 


We are now in position to discuss the expression ( [4.27| ) for the correlation 
C^s of two conserved observables and Qs when some invariance (here the 
A^-invariance) is broken. Depending on the nature of the observables, we 
can distinguish three cases. 


(i) If both and Qs are such that 


S:|A/',K„] = ([Q,Af|>=0 , 


(4.35) 


in particular if and Qs commute with N (as well as with K), their 
expectation values (|4.9| ) are well defined, although Hq is not. Indeed, 
all the density matrices U7Zol7~^ for which Fq{7Z} is minimum, will 
provide the same value for (Q^). The condition ([4.35|) ensures that 


the equation (|4.24|) for d67l/d^s\^=o has a finite solution, that we can 
denote —/d“^Fo ^ : Qs. Due to ( |4.34|) this solution is not unique ; it 
is only defined within the addition of ie [A/*, T^o], where e is arbitrary. 
The condition (|4.35|) , however, implies that this addition is irrelevant. 
Hence the correlation ( \4.2Ti ) is finite and well-defined, in spite of the 
vanishing eigenvalue (|4.34|) of Fq. This conclusion holds in particular 


for the fluctuation AN obtained when Q7 = Qs = N. 


(ii) If Qs, but not Q^, satisfies (U-33i) , the expectation value {Qfi de¬ 
pends on the specific solution Hq of Eq. ( [4.2| ). This feature is character¬ 
istic of a situation with broken iV-invariance. Moreover, the additive 
arbitrary contribution ie[Af, T^o] to d6Il/d^s\^=o also contributes to 
( [4.27|) , so that the correlation between Q 7 and Qs is finite, but ill- 
defined. 


(iii) Finally, if both Q.y and Qs violate the equation (|4.24|) has 

no hnite solution. The fluctuation of Q^, in particular, is infinite, 
consistently with the fact that (O 7 ) is ill-defined. 
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These conclusions hold without any changes for the breaking of invariances 
other than N. 


4.3 Fluctuations and Correlations of Non-Conserved 
Observables 

In this Section we consider the general case of single-quasi-particle observ¬ 
ables which do not commute with the operator K. 


4.3.1 Kubo Correlations 


When the observables do not commute with K, the expression ([4 .161) is 
still the variational approximation for InTr exp[—(3K — ^^ 7 ]- expan¬ 

sion in powers of the ^.y’s, this expression generates cumulants of the Kubo 
type. In particular, while its zeroth and first-order contributions coincide 
with ( |4.7|) and ([4.9| ), respectively, the second-order contribution. 


_ 

^7(5 — ^<57 = 


Tr e 


-pK 


du e 


uK 


Q 'y C 


-uK 


Qs 


/3Tr e-/^^ 


(Q7(4) , (4.36) 


yields the Kubo correlations between the QJs, which are thus given varia¬ 


tionally by ('[4.27I) , that is 




■ ^ 0 


: Q5 


(4.37) 


4.3.2 Standard Correlations 

In order to obtain our variational approximation for the true correlations 
C-ys = ^{Q'yQs + QsQ'y) — {Q'y) {Qs) between and Qs, we have to return to 
the coupled equations of Sect. 3.3, facing their u-dependence and expanding 
them in powers of the sources Sect. 4.1 gave the zeroth-order solution, 
which happened to furnish both the thermodynamical potential and the ex¬ 
pectation values. Fluctuations and correlations are provided by the next 
order. 


We start from the fundamental relation ( 2.20 ) that involves the stationary 
operator 'D{P) and the operator A(^) dehned in (|1.6|) ; more explicitly it 
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reads : 


"'^^sQs -i^-v)'^^sQs 


dip 


TrVm / dwe ^ Q o 
Jo 


-'^^sQs 
Tr V{p) e <5 


(4.38) 


As discussed in Sect. 4.1.2, Eq. ([4.38|) taken at ^ = 0 yields {Q^) as the expec¬ 
tation value (|4.9|) of in the state Voi/d) already determined in Sect. 4.1.1. 


The correlations C^s and the fluctuations AQ^ = are obtained by ex¬ 


panding ( |4.38|) up to hrst order in the ^’s which appear both explicitly and 
through the expansion 

VifJ) ^ Vo{f3) + A{f3) + ...^ X)o(/3) + + ... (4.39) 

<5 

of the stationary state For shorthand we denote with an index 6 the 

coefficient 

dV{P) 




d^5 


(4.40) 


?=o 


of in the hrst-order contribution to Note that, thanks to the vari¬ 

ational nature of the method, only this hrst-order contribution is needed to 
evaluate the wanted second derivatives of the characteristic function. 


By using the formalism of Sects. 3.1 and 3.2, we can rewrite ([4.38|) as 


^ =-q^- |tr2„ (/?) 1^' du Q., 


= -q^ - |tr2„7^[‘^“^(^) f du 

'J 0 








(4.41) 


The correlation C^s is the derivative of this expression with respect to ^s, 
evaluated at .^ = 0. Contributions arise here from the dependences in of 
both 7^b'^l(/d) (or 7^['^“1(/?)) and = J^'y^-yQ-y- We expand as in (|4.39|) the 
matrices 7^W(.^) and TW(,^) around = 0) introducing the same notations 
= dTZ^'^yd^^\^=o and = (9TW/5 ^.^|^=o as in ( 4.40|) (the symbol [i] 
stands for [a], [d], [ad] or [da]). We obtain 


C.yS = 


-^tr2n Q-y'JZ^s'^\p) + itr2n [ Qs , Q-y]'JZo 
-ltT2n (/9) - itr2n [ 

-itr2.Q,{7^'“"^(/?)+7^f'(/?)} . 


(4.42) 
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We now work out the coupled equations of Sect. 3.3 so as to find an 
explicit expression for (|4.42|) . First we note, from the boundary conditions 
( |3.44|) , that we have 


Tt\(3) = 1 , 




\(3) = -Q6 


and hence, using (|3.24D and (117^) : 




ad] I 


S (/3)= 




,[d] 


(4.43) 


(4.44) 


The insertion of (4.44|) into the expression (4.42|) of C^s yields 


C75 — |tr2n — 7^0) + 67(1 ~ 7^0)2.57^0] ~ 

(4.45) 

The hrst term in the right hand side is what would emerge from a naive calcu¬ 
lation using Wick’s theorem with respect to the state Vq{I3). The departure 
from this result, given by the last term, arises from the variational nature 
of the calculation which modihes 'i7(/7) and hence the associated contraction 
matrix when sources are present. 

Rather than evaluating 7i}^\l3) we shall transform Eq. ([4.42|) so as to 
bring the “time” u down to 0 and thus hnd a more explicit expression for 
the correlation To this aim we expand the equations of motion ( |3.52| ) 
for and around the lowest-order contribution (14.3). This leads to 


d7^ 


ad] 


= ^Ko : 


d7^ 


[da\ 


1 . 

2 ' 


= -tKo : 7^ 


dti 2 ’ du 

where the kernel Kq is dehned, for 7^ = T^q + ^7?., by 

Ko:57^ =6[n{n},n] 


[da] 

S 


(4.46) 


= [7^{7^o},57^] + [^ 


dn{TZ} 


, dTZnn' 




(4.47) 


7^=7^o 


The {2n x 2n) x (2n x 2n) tensor Kq is nothing but the usual RPA tensor 
associated with the HFB Hamiltonian 7-f{77} and customarily introduced to 
describe small deviations around the equilibrium state IZq. As in Eq. ([4.22|) , 
the colon symbol in Kq : 5TZ stands for half a trace (with a twist in the 
indices) involving the last pair of indices of Kq and the pair of indices of 577; 
equivalently, Kq : 577 can be considered as a scalar product in the 2n x 2n 
vector space of 577. It is shown in the Appendix (Sect. A.3) that Kq is related 
simply to the matrix Fq which, according to (|4.2CI|) , characterizes the second- 
order variations of the thermodynamic potential Fg{ 77} around TZq. More 
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precisely if we define, for any matrix A4, the tensor Co as the commutator 
with the HFB density matrix TZq, 


Co:Af = , 

the matrix Kq is (in the Liouville space) the product 

Kq = —CqFo = —Co : Fo , 


(4.48) 


(4.49) 


again dehned as in ([4.22| ). The occurence of the simple relation ([4.49| ) and 
its meaning are explained in Sects. 3 and 5 of Appendix A. 

The boundary condition T['^1(0) = 1 implies 


_ ^[da] (0) = 3 ; . 


(4.50) 


A formal solution of Eqs. ([4.46|) , together with the boundary condition ([4.50|) , 
gives 

.huKr 


5 ^ : 3^ , (4.51) 


= e2 “^o ■ 


:3^, 

where y is still unknown. The correlation (^4.42| ) can now be written as 

= —Q^ ■. cosh i/3Ko : 3^ . (4-52) 

With the aim to hnd 3^, using (^4.44| ) and again ([4.51| ), we can obtain 

= [7^o, Q 5 ] = 2 sinhi/?Ko : 3^ . (4.53) 

If Ko had no vanishing eigenvalue, we could eliminate 3^ between ( [4.52| ) and 
(|4.53|) and write: 


= -\Q^ : coth^^Ko : [T^o, Qs] 


(4.54) 


Unfortunately, due both to the vanishing eigenvalues of Co associated with 
all one-quasi-particle operators that commute with ^cx^I-LqCx and to the pos¬ 
sible vanishing eigenvalues of Fq associated with broken invariances (see 
Sect.4.2.3), the operator Kq = —CqFq has no inverse. 

In order to determine the unknown quantity 3^ = 7^^“^(0) entering ([4.52|) , 
we return to the equation (|3.43|) which governs the evolution of Tg{u). 
Eq.( [4.43D gives the boundary condition on while, in agreement with 

Eq. (p.51|) , ^^^“'(O) is related to 3^ through 


ri“^( 0 ) = (l-7^o)-'3^(l-7^o) 


-1 


(4.55) 
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We now relate to by expanding Eq. (|3.43|) to first order 

in the sonrces around the zeroth-order solution (|4.1| - |4.4|) . This yields the 
differential equation 


du 


- = (Eo : 


where the n-dependence has been omitted in and TZ 

the expansion 1-L{TZq + 51Z} ~ TYo + Eq : 51Z, entailed by the dehnition ([4 .191) 


,[da\ 


(4.56) 
We used 


of the matrix Eq. The quantities TZf and 7^^ are given by Eq. (|4.51|) in 


terms of 3^ and as functions of m . In order to transform the left hand side 
into a derivative, we multiply Eq. (|4.56|) by exp |(/3 — u)Ho on the left and 
on the right. As a result, terms of the form e'“^oA7e“’'^° appear on the r.h.s. 
of this equation. Using the fact (see Appendix, Sect. A.3) that the matrix 
So Co in the Liouville space describes the commutation with I3Hq, we rewrite 
these terms as 




SoCo 


M 


(4.57) 


Relying now on ( [4.51| ) we note that, after this transformation, the two terms 
on the r.h.s. of Eq. (|4.56|) depend on u only through exponentials of SoCo and 
Ko. In addition, these terms can be explicitly integrated over u by means of 
the identity 


-SoCo 

d 

du ‘ 


En 


.g±iuKo _ 

-i/3-'iSoCo ±itjKo I r^-1/ 


=F2 —- l)C7^e^5“Ko Q-i('g±|«Ko _ g±i/3Ko)] 


(4.58) 

which is checked by using Eqs. ( [4.21|) and ([4.49|) . Although Co has no inverse, 
the r.h.s. is nevertheless well dehned because Cq ^ is always multiplied either 
on the left side by SqCo or on the right side by CoFq when the exponentials 
are expanded. This amounts to extending to operators the natural convention 
(e“® — = a — b for x = 0. Integrating the transformed differential 

equation for t}°'\u), with the boundary condition = —Qs given in 

(TO) , we End 




(1 


i^SoCoN(^-l 


e2 


mKo 


(1 


e-5^SoCo^(^-lg-i«Ko 


(4.59) 


+2Co ^ (sinh 2/7Ko — sinh ^mKo) ; 3^ 


For M = 0, Eqs. (^4.55|) and (|4.59|) provide the required equation for 3^. We 
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can simplify the result by rewriting (|4.55|) as 


iIpt-Co = 4cosli[|/37io] 3^ cosh[|/37io] 


= -2(sinhiSoCo)C 


(4.60) 


an equality that is derived by using the expressions ( |A.11|) and (|A.19D of Sq 
and Co in the basis where Hq is diagonal. Again the r.h.s. of ([4.60| ) is well 
dehned, despite the vanishing eigenvalues of Cq. We eventually hnd 


= 2Co 1 sinh(i/?Ko) ; 3^ 


(4.61) 


Assuming that Fq ^ exists, Eq. (|4.61|) can be inverted in the Liouville space; 
using the relation (4.49), this gives 






2 sinh |/3Ko 

Hence, from (14.52|) we hnd 

C^s = : (KoCothi4Ko)Fo' : 


(4.62) 


(4.63) 


This hnal result is our variational approximation for the correlations 
or for the fluctuations AQ^. The symmetry C^s = Cs'y is a simple conse¬ 
quence of the relation Kq = —CqFo, whereas this symmetry was not obvious 
from our starting equations ([4.42|) , ( [4.45| ) or ( [4.54|) . The formula ([4.63|) gives 


a meaning to the formal expression (|4.54|) which, due to the vanishing eigen¬ 
values of the commutator Cq in Kq = — CqFo (Eq. (|4.49|) ), was not dehned. 
As for the vanishing eigenvalues of Fq, our discussion of Sect. 4.2.3 still holds 
for (|4.63) as well as for (|4.27|) , and all its conclusions remain valid when 
and Qs are arbitrary single-quasi-particle operators. 

Finally, when one of the observables, say Q^, is a conserved quantity (i.e.. 


[Q^ , A'] = 0), we hnd from (^4.34|) that 


Q 7 : Kq = —Q 7 : CqFo — 0 , (4.64) 

even if the Q.y-invariance is broken (if it is not, we have : Cq = 0 ). 
Then the left-action of the operator KoCoth^/dKo on amounts to a 
multiplication by the c-number 2//?, so that from Eq. ([4.63|) we recover the 
result (|4.27|) of Sect. 4.2.2 for the correlations when one of the two observables 
Q-y or Qs is conserved. 
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4.3.3 Diagrammatic Interpretation 


The expression ([4.63|) can alternatively be obtained by means of a perturba¬ 
tive expansion, where we take as the unperturbed state and \ol^1-LqOl 

as the unperturbed Hamiltonian which generates the quasi-particle or quasi¬ 
hole propagator. Indeed, this expression coincides with the sum of the con¬ 
tributions of all the bubble (or open ring, or chain) diagrams. These dia¬ 
grams, which start from Qs and end up at Q^, iterate the interaction V (see 
Eqs. ( p.27| and (|3.29|) ) with a pair of quasi-particle or quasi-hole lines. A 


partial result [^| was already obtained in the special case when (i) there 
is no broken invariance and (ii) the observables Q.y and Qs are conserved. 
Moreover, the derivation disregarded the difficulty caused by the vanishing 
eigenvalues of the kernel Kq. 

Let us sketch here a proof which takes into account this difficulty and 
applies to systems with pairing. A single line in a diagram carries two indices 
A and A', each of which is associated with either a creation or an annihilation 
operator. Like the contraction matrix TZq, the associated propagator is thus 
a matrix in the 2n-dimensional space introduced in Sect. 3.1 ; it depends 
moreover on the difference in “times” u, with < u < (3. In bubble diagrams, 
the propagator which connects two successive interactions involves a pair 
of lines. It appears as a (2n x 2n) x (2n x 2n) tensor, depending again 
on a single time-difference and having the period (3. As usual, the Fourier 
representation replaces this time difference by an index uj which takes the 
values UJ = 2Tiij3~^m, where m is an integer. Again, it is convenient to 
regard the resulting propagator r(a;) as a matrix in the Liouville space. In 
a basis where TZq and Tio are diagonal, a pair of lines is represented by 




’^OA — 


0)1 


/lom ~ ^OA ~ a; 

The numerator comes from the pair of statistical factors TZq or 1 — T^q asso¬ 
ciated with the lines A and pi ; the denominator accounts for the energies of 
these lines. For = Hqx, which implies TZofj, = TZox, fhe propagator r(a;) 
vanishes, except for a; = 0 in which case it is equal to 


r(AA')(MM')(0) “ 25 a/.j'^/.jA'/37?.oa(1 — T^ox) 


(4.66) 


We can formally rewrite (^4.651) in an arbitrary basis of the Liouville space by 
using the expressions ( |A-26| - [A^) for Sq ^ and the definition ([4.48|) of the 
commutator Cq. This provides 


r(a;) = 
r(0) = 


/3-iCoSo + a; 


(a; 7 ^ 0) 


-/3S 


-1 


(4.67) 

(4.68) 
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The expressions ( |4.66|) and ([4.68|) for the limiting valne r(0) are consistent 
with the valne ( |A.26|) of the metric tensor —Sg ^ for = T^qa- 

Bnbble diagrams iterate a vertex associated with the two-body interaction 
V by means of the propagator r(a;). Regarded as a matrix in the 2n x 2n 
Lionville space, the interaction V can be identihed with the matrix Eg (see 
Eq. (|4.19|) ) of the second derivatives of E{71}, as seen from the expression 
(|3.31|) . Finally, the initial and final vertices associated with the operators 
and Qs (the correlation of which we want to calcnlate) give rise to factors 
and Qs, regarded here as vectors of the Lionville space. Altogether, the 
contribntion of the n-th order bnbble diagram to C^s, evalnated with 'Dq(P) 
as the nnpertnrbed density matrix, eqnals 


C$> = Q, : ;3-‘ |-r(a,)E„]” r(..,) : Qj 


(4.69) 


The snmmation over oj = 2'k\( 5~^ m acconnts for the fact that Q^ and Qs 
are both taken at the “time” n = 0. The ordering of and Qs is relevant 
in the evalnation of the zeroth-order term, where the snmmand behaves as 


— Cg/o; for a; —> ±cx) ; the average over both orderings as in (|1.8|) is achieved 
by taking a principal part in the snmmation over uj. 

The contribntion of all bnbble diagrams to the correlation C^s is thus 




1 


1 -|- r(a;)E| 


-r(a;) : Qs 


(4.70) 


The sum is meant as a principal part for uj 
(|4.67|) of r(a;) we hnd for a; 7^ 0 : 


±cx). By using the expression 


1 -l- r((n)Eg 




/3 ^CgSg -|-a; — CgE( 


■Cg = 


a; + K( 


-Cr 


(4.71) 


where we utilized the relations (|4.21|) and (|4.49|) which introduces Kg. The 
summation over uj = 27ri/9“^ m is then performed by means of 


1 ^/ 1 




(4.72) 


where is meant as a principal part for uj —> ±cx), the value a; = 0 being 
excluded. The expression ([4.72|) approaches a hnite limit for z ^ 0, which 


we take as the value for z = 0. Using now the expression ([4.68|) for r(0), we 
obtain 

C!^5 = Qy ■■ ([-I COth i/?Kg + (/?Kg)-l]Cg + [/?Fg]-') I • (4.73) 
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The first term is defined even if Kq has vanishing eigenvalnes, since these give 
no contribntion to the bracket. The second term may diverge for vanishing 
eigenvalnes of Fq, in agreement with the discnssion of Sect. 4.2.3. Finally, 
nsing again Eq. ( |4.49| ), we hnd 

= Q-y-.ik coth i/?Ko)KoFo ' : Q, , (4.74) 


which is the same expression as onr variational resnlt ( [4.63|) . 

Likewise, if we wish to evalnate the contribntion of the bnbble diagrams 
to the Knbo correlation (|4.36|) , we have to associate with the hnal operator 
Qy of the chain of bnbbles a “time” u rnnning from 0 to /?, instead of fixing 
this time at 0. This amonnts to keeping only the term a; = 0 in the snm 
( |4.69| ) associated with the n-th-order diagram as well as in the snm (|4.70|) 
for the overall set of diagrams. The result reduces therefore to the last term 
of ( [4.73|) , which is identical to our variational approximation ( [4.37| ). 

Thus, both for the standard and the Knbo correlations of single-quasi- 
particle observables, our variational approach leads to the same result as 
the summation of bubble diagrams in a perturbation theory where the HFB 
state is taken as the unperturbed state. The occurence of the RPA kernel 
is consistent with the summation of this set of bubble diagrams, and the 
outcome (|4.63|) shows how the RPA appears naturally in a variational frame 
suited to the evaluation of correlations. The results of Sect. 4.1 show also 
that the bubble diagrams, in our variational treatment, contribute neither to 
the thermodynamic quantities nor to the expectation values (Q^), which are 
both given by the standard HFB mean-field answers. (Diagrammatically this 
corresponds to the summation of the instantaneous part of the self-energy 
insertions.) 

Another variational principle, based on the maximization of the entropy 
rather than on the Bloch equation ( p.3| ) for the characterization of the state, 
has been previously used to evaluate the correlations C^s- While the 


formula obtained for the Knbo correlations was the same as ( |4.37|) , the result 
for ordinary correlations differed from (|4.63|) . Its perturbative interpretation 
also involved the set of bubble diagrams. However, starting from second 
order, only the a; = 0 contribution to each diagram was included (see Sect.4.1 
of [^). In contrast, our result ( |4.63| ) includes the summation over all values 
a; = 27r i j3~^m. This appears as an advantage of the present variational 
principle based on Bloch’s equation. 
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5 


Projected Extension of the Thermal HFB 
Approximation 

5.1 Generalities and Notations 


In Sections 3 and 4 we have dealt with equilibriums in a grand-canonical 
ensemble. However, as discussed in the Introduction, many problems involve 
data of type I that prohibit fluctuations for the associated conserved opera¬ 
tors. For instance, in the case of the canonical ensemble, the state is required 
to have a well-defined particle number Nq. Trial states of the type (|3.4|) that 
were used in Sects. 3 and 4, namely exponentials T of single-quasi-particle 
operators, violate this requirement. We want, nonetheless, to keep the con¬ 
venience associated with the operators T and the Fock space. To resolve 
this dilemma, in conjunction with the variational functional O we shall 
use density operators V involving projections, such as in Eq. ( |1.11|) . For 
instance, a natural choice for a trial state with a well-defined particle number 
is 

V = Pno , (5.1) 

where Pnq is the projection 


Ao = ^ r de 

ZTT Jo 


(5.2) 


onto the subspace with the particle number Nq. As in the preceding Sections, 
the variational parameters associated with the choice (|5.1|) are contained in 
the operator which belongs to the class (|3.4|) . 

We shall consider successively the two types of situations that were ana¬ 
lyzed in the Introduction, (i) The first type (Sects. 5.2 and 5.3) corresponds 
to the case where the choice ( b.l|) involves an operator which breaks a 
symmetry of the Hamiltonian. An example is that of a system in canonical 
equilibrium when pairing sources are added. Then, the operator does 
not commute with Pjvo; although the exact state D = e“^^ does. In this 
case, projections are necessary on both sides of (ii) In the second type 
(Sect. 5.4), no symmetry is broken by the operator This situation is 
illustrated by the thermal Hartree-Fock approximation. Then, the operator 
does not violate the A^-invariance and V can be defined with only one 
projection on either side of this projection is needed because acts 
in the unrestricted Fock space. 

Whatever the case, projected trial states are not easy to handle in the 
form (^). In particular, they do not satisfy directly Wick’s theorem which 
has been an important tool in Sects. 3 and 4. In order to retain the simple 
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properties of the exponential operators T, we wish to express each projec¬ 
tion as a sum of such exponentials. This is feasible when the projection is 
associated with some underlying group of operators having single-quasi- 
particle operators as generators. In this case the projection can be expressed 
as a weighted sum of the operators over the elements g of this group. 
For instance, the projection in ( b.2|) is built from the group elements 
j^ig] = ^ 0 '^ which are indexed by the parameter 0 < 0 < 27r. As 

a consequence, our projected trial states will appear as sums of simple expo¬ 
nential operators of the T type, each of which is easy to handle. 

Likewise, if we only wish to eliminate the odd (or even) particle-number 
components, we should introduce a projection on even [rj = 1 ) or odd {rj = 
— 1 ) numbers, which can be expressed as the sum 

P, = \{l + Ve^-^) , v = ±l (5.3) 


over the two elements of the group e^®^ with 9 = 0 and 9 = n. 

The projection on a given value M of the ^-component of the angular 
momentum, 


Pm = 


1 

2 tt Jo 


( 5 . 4 ) 


is similar to ( |5.2|) . To project in addition on a given value of the total angular 
momentum J, we should introduce the rotation operators = R{il) in 
the Fock space ; here the group index g denotes the three Euler angles hi, 
and R{Vt) has again the form of an exponential of single-particle operators. 
Integrating over the Euler angles with the rotation matrices D'ljp^{Q) yields 
the operator 


2 7 -I- 1 /■ 

PMK = ^frJ'iaDi;ja)R(n) = J2hJM}{',jK\ , ( 5 . 5 ) 


where I 7 J M) denotes a basis in the Fock space labelled by J, M and another 
set of quantum numbers 7 . When K = M, Eq. (^.5|) dehnes the projection 
on given values of J and M. 

The projections P± over the spaces with a given spatial parity can be 
expressed as 

P± = i , (5.6) 

where the single-particle operator 


N_ 


I ^ y d^r [fjlir) - fjii-r)] [V'<,(r) - 


(5.7) 
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counts the number of particles lying in the single-particle states which are 
odd in the exchange r, a —> — r, a (the factor | in Eq. (|5.7|) accounts both for 
the normalization of the basic antisymmetric states [t/’J(r) — ' 0 l(—r)]| 0 )/A /2 
and for the double-counting in the integration over r). 

In order to encompass formally all these situations, we shall denote by 


P= 


(5.8) 


the decomposition of our relevant projection as a linear combination of the 
group operators The symbol /^ • • • denotes a weighted sum or inte¬ 

gral over the elements of the group. In particular, for ( |5.2| ), it stands for 
(27r)“^ ... while for (^), it stands for \ I]j=o,i ■ The operators 

have the form (^.41) of an exponential of a single-quasi-particle operator, 
which we shall denote as — The coefficients and £1^1 de¬ 

pend on the group index g ; for instance, in (p.2|) , we have = —i9n/2 and 
£[^1 = —i9Af {0 < 9 < 2 tt), and in (|0|), = —ij 11 x 1/2 and = —iinAf 

{j = 0, 1 ), where A/” is the 2n x 2n matrix dehned by ([4.12|) . 

As in Sect. 3.1, we denote the 2 n x 2 n matrix exp(—£[®]) as £[^1 and 
the scalar exp(—as We shall have to deal with products of two to 
four T-operators, such as . from the algebra 

properties (|3.23| - |3.2^) , these products have still the same T-form and they 
are characterized by matrices such as 


q-lgiig'] ^ 'j-[g]'j-[d.]'j-[g'] 


(5.9) 


and scalars = p 9 ]j-ld\^[g']_ There is a contraction matrix 'R/^dg'] associ¬ 
ated with the matrix T^sdg'] through ( |3.8|) ; as T^^dg']^ depends on the 

ordering of the original operators. 

The equations that we shall derive in the next Subsections can also han¬ 
dle cases which are more general than pure projections. For instance, in the 
description of rotational bands of axial nuclei, it is often a good approxima¬ 
tion to assume that the quantum number K associated with the component 
of the angular momentum on an internal axis of inertia has a well-dehned 
value Kq. Then rotational states with given J and M values can be sought 
with the trial form 

V ^ PiPirn , (5.10) 

where the operators Pmk dehned by Eq. o are not projections when 
M 7 ^ Kq. Thermal properties of triaxial nuclei which involve sums over the 
K quantum number, can also be handled by our method. Another example 
is the case where a symmetry (such as N) is broken. Suppose that we still 
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want to work in the grand canonical ensemble bnt also want to snppress the 
off-diagonal contribntions between snbspaces with different values of N. We 
are led to consider operators of the form 


Y^PnDPn = 

N 


n27V 1‘2'K 

1 r^TT 


d(j)' , 


N 


(5.11) 


27r Jo 

which are also amenable to our treatment. 


5.2 The Variational Ansatze and the Associated Func¬ 
tional 

Now we must choose the trial classes for the operators P{u) and A{u) to be 
inserted into the functional ( ^.41) which, at its stationary point, supplies the 
characteristic function of interest. For the operator T>(u) we already made 
the choice 

V{u) = Pf^'^{u)P = [ , (5.12) 

J g g' J gg' 


where T^^{u) is defined by Eq. (|3.16|) . For the operator A{u) we keep the 
same choice, 

A{u) = f^°^{u) , (5.13) 

as in Eq. (|3.17|) . Again, as in Sect. 3.2, the quantities P^{u), /[“^(n), 

£[“!(«) form a set of independent variational parameters. The exact density 
operator D satisfies PD = DP = D, and when defined in the full Fock 
space the operator K commutes with P. Therefore, it makes no difference if 
we project D{u) or A{u) in the functional (p^); the choice D{u) = T[‘^(m), 
A{u) = PT[“1(m)P is equivalent to (|5.12| - |5.13|) . 

The results of Sect. 3.2 are readily adapted through a mere replacement 
of by pbiia ] amj a summation over the group indices g and g'. Thus the 
normalization Z, given by 


Z = Tr A{u)V{u) = / Z39 (m) 

J gg' 


(5.14) 


with 


\nZ99'(u) = -[/[“'(M)+/['^'(M)+/[3'+/[^'l] + itr 2 nhl[l + rW(M)r[®lr['^](M)T[^?'l] ^ 

(5.15) 

replaces Z as dehned by (|3.25|) . We recall that, according to Eq. (|2.17|) , Z 
does not depend on u when A{u) and D{u) are stationary solutions. 
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The action-like functional (p.32|) , is now replaced by 
r , ( 


J = - / du 

Jo 


Zaa 


9 9 


+ -tr 2 n 'Ridg'ag] ^ 1 | ^ ^{7^N'^9']}] \ _ 

2 dxi 2 j 

(5.16) 

The boundary term = Tr A'D{/3) is obtained from Eqs. ( 5.14| - |5T5|) by 
letting u = fJ in and T^‘^{u) and replacing and £[“](«) by the 

given quantities and (see Eq. (|3.14|) ). We have used the analogue 
of the property (|3.33|) to transform ^Ji[ 9 dg'a]>j-[g\ j^Wag]^ Moreover, 
from the commutation [K ^ T^^l] = 0 and in agreement with ( 4.31 ), we find 
that for any TZ the energy (|3.31|) satisfies 


E{T^9]^T^9] = E{n} , 


(5.17) 


a property that we have also used to write instead of 

(similarly, one has 

The expression (^.16|) exhibits the fact that we might as well have pro¬ 
jected A{u) instead of Viu). Indeed, taking into account the property (^.17|) , 
this expression does not change under the replacement of gdg' by d and of a 
by g'ag. Actually, we can readily check that if one projects A instead of V, 
Eq. ( p.32| ) is directly transposed into (|5.16|) through the replacement of a by 
g'ag. 

Likewise, by an obvious extension of ( 3.35 ) where d is replaced by gdg', 
we find for the functional I the alternative form: 


J = - Z{(3) + Z(0) + [^ du [ Z^^' (- 

Jo Jgg' \ 


d/[“l(M) 

dw 


+ -tr2„ T[“1 ^ n^^^adg'] "j . 

( 5 . 18 ) 

The expressions ( |5.16| ) and ( ^.18|) generalize (|3.32|) and ( p.35|) through the 
averaging over g and g' with the weight Z^^ /Z. They lead us to introduce 
the matrices 

z3 9'giWag]^ ^[a.d.] ^ ^-1 f ^5g'^[agdg'] 

'Jgg' 'J 0 q' 


gg 


and the pseudo-energies 
2 (da) ^ f E{71^‘^^'°'^^ 


E(ad) ^ J” 


9 9 


-1 /■ ^9 9'^{7^k9rf9']| ^ 

J 9 9' 

(5.20) 
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(remember that and that E{Tl^°'^'^^'^} = E{TZ^3'°'^'^^}). 

With these notations we can rewrite (|5.16|) and (|5.18|) in a more condensed 
fashion, formally closer to the fnnctionals ( ^.321 ) and (|3.35| ) : 

d/[‘^ 


I = 


dn Z 




du 

’ 1 ^ _l_ l|'_g'W _g'Mj 


dn 


r/3 


d/[°l 

dn 


J = - Z{p) + Z(0) + [ du Z 

Jo 

HTW 1 _ _ ^ 

/V' - 


1 


+ -tr 2 nT[“ 
2 


(5.21) 


dn 


where '^['^•“•1, '^[“•'^•1, and E^°‘^'> replace 7^1'^“!, 7^1“'^!, E{TZ^‘^°‘^ and E{TZ^°‘^} 
while Z replaces Z. An important change between the present fnnctional and 
the one of Sect. 3.2 lies in that differs in general from E{'JZ^‘^'°"^ (as well 
as from E{'R}'^'°‘^) when the operator K contains a two-body interaction. It 
implies that the variational ontcome Z^-^^ yielded by the fnnctional (|5.21|) will 
differ from that which would be obtained by projecting the optimum state 
found in Sect. 3 ; this holds even for the case {A = 1 ) of thermodynamic 
potentials. 


5.3 The Coupled Equations 

As in Sect. 3.3 we find the equations governing the evolution of /[‘^(m) and 
(m) by requiring that the functional (|5.16|) be stationary with respect to 
/[“](m) and T[“1 (m). The resulting differential equation for /['^I(m), 


^ = ^tT2n ^ 


(5.22) 


is similar to Eq. (|3.36|) , with '^['^•“•1, and replacing 7^1'^“!, 

and E{7^[“'^1}. The differential equation for T['^1 (m) differs from Eq. (|3.37] ) by 

an explicit summation over the group indexes. It reads : 


f dTt'^l 


J ^9 9''J-[9](^l _ I_ 


_|_ \^Wa9\'J-[d\ _|_ > (1 — = 0 


(5.23) 


with 


JdA =7^{7^W} + 


ky 


1 - 7^W 


(5.24) 
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where i stands for dg'ag or g'agd and where the c-nnmber ^ is dehned by 

kd^ = ^ jtran ] 

(5.25) 

To derive the equation (|5.23|) we made use of (|3.51|) , of ( b.22| ) and of the fact 
that 7i{7?.} satishes for any IZ the relation 


n{T^ 3 ]gir^ 9 ] 1 } = T'^ 9 \n{n}T^^^ ^ , ( 5 . 26 ) 

a consequence of the dehnition (|3.38| ) and of ( b.l7[ ). Note that Eq. ( |5.22| ) 
does not imply the vanishing of the integrand of (|5.16|) before summation 
over g and g'. 

Likewise the differential equation for /[“!(«) is obtained by varying the 
functional (|5.18|) with respect to : 


1-1 




1 _ htN 1 


( 5 . 27 ) 


while the equation for is provided by varying (|5.18|) with respect to 

(m) for 0 < M < /? : 


9 9 ‘ 


Z3 9 'q-W] _ qila9dg']'j 
1 


■ dr [“1 

du 


(5.28) 


^'jqla9d9']q-[a] _|_ q-[a]^lgdg'a]-^ m1 “ Jlld^g'a]-jq-[g] _ g 


where 'H\°-9<^g'] and 'H'^sdg'a] are defined according to Eq. (|5.24D with i standing 
for agdg' or gdg'a and with the c-number ^ replaced by k^9' given by 

if _i _ rlT[“l 

k9/ = - I -tr2n [T[“] (n^'^9d9'] _ g^ia.d.] ^ __ ] 

_i_ 

(5.29) 

Note that disappears from Eqs. (|5.23| ) and ( b.28|) if we divide them 

by Z. These equations thus couple only the matrices and Once 

they are solved, and are obtained by integration of Eqs. (|5.22|) and 
( 5.27|) . In Appendix B, using the Liouville formulation and the notation 
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( |5.19| - |5.20|) , Eqs. ( |5.23|) and (|5.28|) are rewritten in the forms ( |B.10[ - [B?TT|) 
which look similar to Eqs. (|3.37|) and ( p.43|) . 

Finally, in agreement with p.l8| ), the stationarity of ( |5.18|) with respect 
to and T['^l(/?) provides the equations 

Z{f3) = Z^^^, , (5.30) 

which are satished by the same boundary conditions as in Sect. 3 : 

. (5.31) 

The boundary conditions at u = 0 are /['^^(O) = 0, T[‘^1(0) = 1 as in Sect. 3. 
Again, all our variational quantities (/[“^(m), /['^^(m), T[“1(m), T^'^{u)) depend 
upon the sources through the boundary conditions (|5.31|) . 

Despite their somewhat complex form, the coupled equations (|5.23|) and 
(|5.28|) entail some simple consequences. First, using ( b.22|) and ( b.27|) , one 
notes that Z^^ varies with u according to 

d 


dn 


-\nZ39' =k{9_^9J!' 

Id 


_c[u) 


(5.32) 


where the function c{u) is given by the average 

-c{u) = -X [ Z9 9' ^tr 2 „ hi(l - 

2Z J gg' du 

Summation over g and g' of Eq. (|5.32|) leads to 


(5.33) 


(5.34) 


which was expected from (|2.17| ). Thus the ultimate quantity of interest, the 
characteristic function ip provided by the stationary value of InX, is equal to 
\i\Z{u) for any u. One can also check from ( |5.23|) , ( b.28| ) and (|5.32| ) that 

= 0 , (5.35) 

du 


which is a special case of the general relation (p.l5|) . The relations (|5.34|) 
and (|5.35|) generalize (|3.48|) . 

When the operator A is hermitian, the coupled equations (^.23|) and (|5.28|) 
admit solutions with hermitian matrices and T[“1 (m), and the result¬ 

ing c-numbers l^‘^{u) and /[“!(«) are real. This is seen by grouping in all 
the equations the terms corresponding to respectively with the 

associated terms corresponding to jg ghown at the end of 

Appendix B. 
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5.4 Unbroken P-Invariance 

The above approach is general and covers situations in which the trial oper¬ 
ators T^^{u) and do not commute with the projection P. It applies 

for instance to hnite systems of fermions in canonical equilibrium in which 
pairing correlations are operative. In this case the operator K denotes the 
Hamiltonian H alone without the term —piV, the projection P is given by 
( |5.2| ) and the independent quasi-particle trial state T^^{u) as well as the trial 
operator do not commute with N, nor with the operations = e'^^ 

of the associated group. 


5.4.1 Notation and General Formalism 

As discussed in Sect. 5.1 and in the Introduction, there are other situations 
which require projection, although the P-invariance is not broken by the 
variational approximation. In these cases and commute with 

all the group elements and therefore with P. Hence only one projection 
is needed on either side of T^^{u) and, instead of ( b.l2|) , the trial state can 
be written as 

V{u) = P (u) P = f {u)P = P f (u) = [si f ['^ (u) . (5.36) 

A single summation over the group index g is sufficient in Eq. (|5.14|) : 

Z = Tt A{u)V{u) = [ Z3{u) , (5.37) 

while Eq. (|5.15|) ) becomes 

InZS(M) = -[/[“l(M) + /['^l(M) + /[sl] + itr2„ln[l + r[“l(M)r[s]r['^l(M)] . (5.38) 


Likewise, the double averaging in (|5.19|) or (|5.20|) is replaced by a single one 
over g with the weight Z^jZ \ for instance, we have : 


= ^-1 f ^ = ^-1 f ^9 7^Ms] . (5.39) 

•Jg Jg 

The expressions of the functional X and of the coupled equations are 
readily derived from those obtained in Sects. 5.2 and 5.3 by noting that the 
matrix Plsl commutes with pl'^l [u) and pl^l {u) and that the factor pls 1 can 
everywhere be replaced by pl^l = 1. Since g' disappears and the ordering 
of g with respect to a and d is irrelevant, there remain only two (instead of 
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four) different contraction matrices: = Jiigda] t^Ms] _ 

j^iagd] _ j^[gad] _ functional ( ^.161) then simplifies to 


J = 


dn 


d/['^'( 


u) 




du 

dT^'^l 1 


(5.40) 

As before, the stationary value of InX (our variational approximation for the 
characteristic function) is equal to InZ which does not depend on u. 

Accordingly, the coupled equations of Sect. 5.3 take a simpler form. Note, 
however, that because in general (Eq. (^.25|) ) differs from (Eq. (|5.29|) ), 
the matrix T-C^sdOa] jg different from (different from 

(Eq. (F^) ). 

By combining the forms taken by Eqs. (p.23|) , ( b.28| ) and (|5.32|) in this 
commutative case and by using (|3.51|) , we hnd that the averaged contraction 
matrices (|5.39|) satisfys 


du 2Z . 

^^{da) _ _^ 

du 2Z 


? 


(5.41) 


This was expected from the general equations (|2.32|) and (|2.33|) , since their 
validity conditions are satished when the operators T[‘^1 (u) and T[“1 (u) com¬ 
mute with P. While Eqs. ( p.52|) , which are appropriate to the unprojected 
variational treatment, had the same form as the time-dependent HFB equa¬ 
tion (with an imaginary time), Eqs. ( 5.41 ) deal with the larger set of matrices 
j^ldag] jiladg] gg^ involve a Summation over the group index g. The ap¬ 
parent deconpling of Eqs. (^.41|) is again hctitious because the boundary 
conditions cannot be expressed explicitly in terms of and . 


5.4.2 Partition Function, Entropy 

In the rest of this Subsection, still assuming that the invariance is not 
broken, we focus on the case A = 1 for which our variational principle pro¬ 
vides an approximation for the exact partition function 

Zp = Tr Fe“^^ . (5.42) 

We have shown in Sect. 2.3 that by shifting appropriately the operator K, 
the variational principle ( |2.4| ) used with A = 1 also yields the characteris¬ 
tic function and therefore the cumulants of the one-body operators which 
commute both with K and with the projection P. 
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As seen in Sect. 2.4, the boundary condition .4(0) = 1 allows us to restrict 
the trial classes (|5.12|) and (|5.13|) to 


Vniu) = Pe 


-uHn 




(5.43) 


the subscript 0 refers to the zero values of the sources The trial operator 
Hq is a (w-independent) single-quasi-particle operator commuting with P, 


Hq — hf) -\- 'H-o cx. , 


(5.44) 


which depends on the parameters Hq, a real number, and TYq, a 2n x 2 n 
hermitian matrix obeying ( |3.41|) and commuting with The choice (|5.43|) 
entails 

Vo{u)Ao{u) = Pe-/^^o = Ao{u)Vo{u) = MP) , (5-45) 


and the w-integration becomes trivial in the functional (p.4|) , which reduces 
to 

J = Tr pK) 


= 


1 + + \(3tl2n HoK - 


with 


InZo^ = -pho - + itr2„ ln[l + , 


(5.46) 


(5.47) 


where the trial quantities and Tio must be determined variationally. In 
the generalized contraction matrices 




j’la] 




e/3Ho + r[5l l + Po(^[^'-l) ’ 


(5.48) 


the matrices TYq, and IZq all commute; for the unity element of the 
group (T[°] = i ; = 0, = 1), the matrix 77 q given by (^.48|) reduces 

to the contraction matrix TZq of the usual finite-temperature HFB formalism 
(Eq. (H)). As in Eq. (|2.43|) , we recognize in the simple expression (|5.46|) 
the standard variational principle associated with the maximization of the 
entropy compatible with the constraint on {K). If one introduces the quan¬ 
tities 




n= z?,nl, E=^ z^EiK] , 


(5.49) 


the functional ( b.46| - |5~47|) takes the more compact form 


I = Z 


1 -|- /3ho -l- ^PtT2n 7fo77 — PE 


(5.50) 
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Within the selected variational space, the best approximation for the 
exact partition function Zp (Eq. (|5.42|) ) is supplied by the maximum of X 
since the r.h.s. of (p.46D is always smaller than Zp. The maximization with 
respect to the number ho yields 


hn — + E 


(5.51) 


This equation ensures that, according to the discussion of Sect. 2.4, the 
maximum of X reduces to Z. Requiring now the functional X to be maximum 
with respect to the matrix Tio provides 


Jo 


with 


where 


w = + 


k9 




= --tT2 n Hoin^o -'E) + E{nf,} - E 


(5.52) 


(5.53) 


(5.54) 


The equation (^.52|) , together with the definitions ( b.53|) and ( b.54| ), deter¬ 
mines the trial matrix TYq to which TZq is related through (^.48|) . The A-basis 
where TYo and IZq are diagonal (with eigenvalues TYoa and respectively) 
is thus dehned by the equations 

(5.55) 

■'3 

while the diagonal eigenvalues of the matrix T^q are given by 




{k» jW{7?|}aa(1 - nU) + nx(E{K) - £)} 

(5.56) 


where M is the symmetric matrix 

AfA„ = /{^a„KSa(1 - KSa) + \Kx(n, - K,)} 


(5.57) 


Therefore the projection generates a more complicated relation (Eqs. (|5.52|) 
and (|5.48|) ) between Tio and the T^q’s than the relation between 7-fo and IZq 
(Eqs. ([4.1|) and ([4.2|) ) given by the grand canonical HFB formalism . 

The equations (|5.52| - ^754D can also be obtained directly from the evolution 
equations ( |5.23| ) and ( b.28| ) when one takes into account the simplifications 
induced by the choice (|5.43H5T4^ and by the commutation properties result¬ 
ing from Eqs. (|5.41|). Indeed, k^ is then the n-independent common value 
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of the quantities introduced in (|5.25|) and (p.29|) , while the 

matrices dehned in ( ^.241 ) also become n-independent and all reduce to 

q^lagdo] _ q^lQagd] _ q^[dOag] _ q^[gdOa] _ 'j^g 

Several consistency properties can be derived from the formulation above. 
When the functional ( 5.501) is stationary with respect to the number ho, the 
approximate entropy S associated with the approximate normalization Z and 
projected HFB energy E by the thermodynamic relation 


S = \xiZ + (3E 


(5.58) 


(a specialisation of Eq. (|2.28|) to the restricted variational spaces dehned by 
Eqs. (|5.43| ) and (|5.44| )) is also related to the approximate density operator 
'Do{{3) = by the von Neumann relations ( p..2| ) and ( ^.45| ). Using 

Eqs. (|5.47|), (|5.49|) and (|5.51|), one obtains the more explicit form 


S = -|tr2„ [7^1n7?.o + (1 - 7^)ln(l - TZq 


+ln/ e-'“[det{l + - l)}]^/^ 


(5.59) 


As expected, S reduces to the HFB entropy S'{7^o} (Eq. (E3» when the 
group of operators is restricted to its unity element. Furthermore, in 
accordance with Eq. (|2.26|) , the stationarity of (|5.46|) with respect to the 
matrix Tio ensures that the equilibrium energy, —d\nZ/dj3, is equal to E. 
According to the discussion of Sect. 2.2, this guarantees that, in agreement 
with thermodynamics, the projected entropy S also satishes the relation 
(12)^) where Eq = -T \nZ. 

Remark : Through the Ansatz (|5.43| ) we have restricted the trial class by 
choosing an exponential u-dependence for the operators Viu) and A{u) and 
by relating their exponents. Therefore, it is not obvious that the maximum 
of X as given by Eq. ( b.46| ) is equal to the stationary value of the functional 
(15.401) which is dehned within the more general class (|5.12| - ^T^ of operators 
whose M-dependences are not limited a priori. In the light of Sect. 2.4, the 
commutation of P with and is sufficient to ensure this equality. In¬ 
deed, if we associate a projection P to A{u) as well as to V{u), the conditions 
i) and ii) of Sect. 2.4 are both satished. 

Nevertheless let us give a more explicit proof along the line already fol¬ 
lowed in Sect. 4.1 for the HFB case. We have to check that the Ansatz 
( |5.43|) with ho and Tio given respectively by ( b.51j) and (|5.52|) satishes the 
set of coupled equations (^.22|) and (|5.23|) , or (|5.27|) and (|5.28|) which express 
the unrestricted stationarity. We hrst recall that, for the choice (|5.43 ), the 
matrix 77['^“^^(n) is u-independent and given by ( |5.48|) . As = uho and 

/[“](«) = {(5 — u)ho, the functional (|5.40|) is identihed with (|5.46|) and both 
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equations (|5.22|) and ( |5.27|) are identical with (p.51|) . The equation (|5.23|) 
then becomes 


Z^qTZI f - ^e“^' 


(1 - 7^g) = 0 


(5.60) 


Moreover, when TIq is solution of Eq. ( b.52|) (which implies Eq. (|5.55|) in the 
basis where Tio is diagonal) one has the commutation relation 


Ho 


z^nl 


ni) 


= 0 


(5.61) 


so that Eq. (|5.6UD is satished for all values of u. This last property can also 
be demonstrated by means of Eqs. (|5.41|) . Likewise, within the replacement 
of M by /? — M, the stationarity condition ( b.28| ) also has the form (|5.6(]|) and 
this is a consequence of (^.52 ). Therefore, in the case A = 1 and unbroken P- 
invariance the result of the optimization of the functional (|5.16|) is identical to 
the approximation Z of the exact parition function Zp (Eq. (^.42|)) obtained 


by the maximization of (|5.46|) . 


66 

































6 Projection on Even or Odd Particle Num¬ 
ber 


In this Section we consider finite systems of fermions, such as atomic nuclei 
or metallic grains, for which the Hamiltonian favors pairing correlations. 
In such cases there exists a qualitative difference between the properties of 
systems with even and with odd numbers of fermions, whereas systems with 
N, N±2, N±4 ... particles display strong similarities. At zero temperature 
the HFB approximation, which is known to provide a good description of 
pairing correlations in nuclei, involves a sum of components with the same 
particle-number parity. On the other hand, at non-zero temperature, any 
number, odd and even, of quasi-particle configurations are allowed by the 
standard HFB solution. As a consequence the HFB thermal equilibrium state 
has both even and odd particle-number components. A completely realistic 
description of finite systems with fixed particle number A^^o would require a 
full projection as defined by Eqs. 


(See Chapter 11 of Ref. ||I0| for 
the litterature, starting with Refs. |P0|| , about the particle-number projection 
at zero temperature.) Nevertheless, the similar properties of systems with 
the same particle-number parity suggest that a first significant step towards 
this goal consists in projecting the independent quasi-particle trial state onto 
the space with even or odd particle (or equivalently quasi-particle) number 
by means of the projection Pjj defined in Eq. (|5.3| ). 

Using the formalism of Sect. 5 we want therefore to approximate varia- 
tionally the characteristic function (PniO of ^ system with a given parity of the 
particle number (or V-parity). In particular, for A = 1, we shall thus obtain 
an approximation for the number-parity projected grand-partition function 






( 6 . 1 ) 


where the prime indicates that the sum runs only over even [r] = 1) or odd 
[t] = —1) values of N. We still enforce the particle-number average by means 
of the chemical potential /i entering the operator K = H — p,N. 

In the specialization of the equations of Sect. 5 to the present case, the 
sum Jg in the general formula ( b.8|) reduces to a sum of two terms: 

= i(fM+77fW) = Ki + r^U-'^) . (6.2) 

In the 2n x 2n representation, we have simply 

/[o] ^ Q ^ q-m ^ I . fi] ^ _ii7rn , rW = = -1 . (6.3) 
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The commutation of and with any matrix T of the algebra intro¬ 
duced in Sect. 3.1 is obvious. This reflects the fact that all operators of the 
single-quasi-particle type ( |3.12|) , or of the trial forms or given by 
(|3.4| ), commute with and and with the projection because they 
involve only an even number of creation or annihilation operators. Thus, 
although a BCS type of state ( p.l6|) breaks the particle number invariance 
e*®^, it does not break the N-parity invariance and we can work with 

the variational spaces 

P(u) = , A{u)=f^°^{u) , (6.4) 

where only one projection is needed. 


6.1 Characteristic Function in the Particle-Number Par¬ 
ity Projected HFB Approximation 

In the present Section we shall therefore rely on the formalism developed 
in Sects. 5.3. The reader who has chosen to skip Sect. 5, is requested to 
accept Eqs. (|5.37| - ^.40D which, after the replacement of by a sum and the 
simplifications associated with Eq. o. furnish the expression of the func¬ 
tional corresponding to the iV-parity projection. The matrices and 

-Ridag] 

(u) which occur in the functional (|5.4CI|) can now be written explicitly: 








2PW 


(6.5) 


with i = ad or da and with or Rida] (^legned as in Sect. 3.2. The weight 
factors Z^, given by Eq. ( 5.38|) , take now the compact form 




- 1/2 


Z^ = rZ\ r = {det2n[A(l - 2p[“‘'l)]}'^' 
and hence the normalization Z^j (Eq. (^.37|) ) reads 

Z, = |Z°(l + 7/r) . 


( 6 . 6 ) 

(6.7) 

( 6 . 8 ) 


According to the general formalism of Sect. 2, the evaluation of the ap¬ 
proximate characteristic function (p{^) = InZ^ = In TrA(,^)'D(/9) requires the 
solution of the coupled equations which determine 'P(m) and A{u). To this 
aim we will use Eqs. ( b.23 - F^2^ , (|5.28| - |5.29|) and (|5.32|^.33|) with the follow¬ 
ing simplihcations : i) suppression of the integration on the group index g' 
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which is everywhere set to 0 = i), (ii) replacement of Jg by 

the sum J2j with j = 0 or 1, (iii) commutation of with and 
whenever it is required. Thanks to these simplihcations we can express the 
evolution equations in a more explicit form than in Sect. 5.3 since the group 
elements are now trivial. For instance, Eq. ( b.23|) for writes: 


^r^d\ I 

—^- rjr 

au 




where we have introduced the generalized HFB Hamiltonians and 


defined by 

^ ^[ 0 a 0 d]_^^ 


-^[Oald] 


27^M _ 1 27^H1 - 1 

nf"'^ = -_ ^[dOai] 1 


( 6 . 10 ) 


2 ']l[da 


1 


27^[da 


in terms of the matrices given in Eq. (|5.24|) . Taking into account that the 
quantities ® defined in (|5.25|) now satisfy the relation 


= -r]rk 


01 


( 6 . 11 ) 


with 


'^d — 


2(1 + rjr) 


2tT2n [Tt"] 


i7^[H(i 7^F«]) dT^d] 


E{-. 


^[da 




1 du 
7^ Ml 

one can also express the matrices and as 


( 6 . 12 ) 


- nm - <1 r ^ H{^} ^ 


27^W - 1^ 27^W - 1 27^W-l 


, (6.13) 


where i stands as above for ad or da. Likewise, the evolution equation ( |5.28| ) 
for Tb](M) can be put in a form similar to (|0|). In Eqs. (|0| - |6.1CI|) and (|6.13|) , 
the subscript d in the matrices accounts for the fact that the c-numbers 
and and therefore the matrices and 7^^*^ are in general different. 
For g = g' = 0, the specialization of Eqs. ( |5.32| ) and (^.33|) to the A^-parity 
projection yields a compact form for the evolution equation of lnZ°, 


dlnZ° gr 1 d7^['^“l, 

du 1 + gr 27^[‘^“1 — 1 du 


(6.14) 
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in 


terms of (or of This equation could also have been obtained 


from the general property dZ^/dw = 0, using Eq. (|6.7|) and the derivative of 
r drawn from (| 6 . 8 D . To complete the set of equations which may be useful 
for the calculation of the characteristic function, we can rewrite Eqs. (|5.41|) 
in terms of the matrices and For the evolution of we get: 


d7^[da 

du 


d7^[da 


rjr 




2n^da] _ I 

= -\[n. 


du 2 n^da] _ I 

(da) 


2n^da] _ I 


du 


(6.15) 

where dhiZ^/du is given by Eq. (|6.14|) . A similar equation, with da replaced 
by ad and a plus sign before the commutator in the right hand side, holds for 
-^[ad] _ 'ppggg equations for TZ^^da] -^[ad] are coupled through the boundary 
conditions, as were Eqs. (|3.52|) which they generalize. 

Remark : Due to the square roots occuring in Eqs. (|0| - |6.7|) , we must 
carefully £x the determinations which control the signs of and Z^. In 
Eq. ( p. 6 |) we can rely on the fact that, for small values of the sources, the 
contraction matrix 'R}°'d^ is nearly real, lying between 0 and 1 , which makes 
the determination unambiguous. In Eq. ( |6.7| ) we have introduced the factor 
M to account for the term —of ( b.38|) , but both the matrices M and 
1 — have positive and negative eigenvalues ; moreover, they do not 

commute when pairing takes place. In order to overcome this difficulty and 
to dehne unambiguously the sign of r, we proceed by continuity, starting 
from the unperturbed state without pairing. In this case, when M and 
are simultaneously diagonalized, the two diagonal blocks of A/'(l — 27^ 
are identical, due to ( 0 ). and provide twice the same factors in (|6.7|). We 
then choose the determination of r as 


r = ^(l- 27^[“‘'lA) , (6.16) 

A 

where are the diagonal elements of lying in the block M = 1. 

Hence, r is positive if the number of eigenvalues of the single-particle Hamil¬ 
tonian that are smaller than the chemical potential /x is even, while r is 
negative if this number is odd. In the low temperature limit, for A = 1 and 
arbitrary /x, one obtains for the grand potential Fq = —(3~^ InZ^ the expected 
result, namely the minimum of {E — fi{N)) for either even N (77 = -|-1) or odd 
N [f] = —1). When the interactions are switched on and pairing takes place, 
we note that within the determinant of Eq. the particle number ma¬ 
trix Af can be replaced by the matrix Af (associated with the quasi-particle 
number) which commutes with We can thus dehne r by continuity. 
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starting from ( |6.16| ) and keeping the same form ( p.l6| ) where the a are 
now the eigenvalues of associated with the subspace M' = +1. 


6.2 The Partition Function in the A/^-Parity Projected 
HFB Approximation 

If one is only interested in the evaluation of the iV-parity projected grand 
partition function (obtained for A = 1 ), or more generally of the char¬ 
acteristic function for conserved single-particle observables (obtained for 
A = with , K] = 0), it is not necessary to solve the 

coupled equations of Sect. 6.1 which are suited to the evaluation of the most 
general characteristic function. 

As in Sect. 5.4.2, we introduce a trial operator Hq of the form ( b.44|) which 
depends on the c-number and the 2 n x 2 n matrix Tio; for the restricted 
variational spaces, we take 


Va{u) = F, , A(ti) = 


(6.17) 


The matrices and are therefore equal to exp(—wT-fo) and 

exp(—(/3 — M)7-fo), respectively, while the w-independent contraction matrices 
and are expressed in terms of Tio by 




1 + e^^o 


— 7?.n 


(6.18) 


We shall also need the specialization of the matrices TZq, dehned in Eq. (^.48|) , 
to the simple projection group associated with : 


TZq — TZq , 


7^1 _ ^0 

° 27^n - 1 


X — 


(6.19) 


As in Sect. 5.4.2, the c-number ho is given by Eq. (|5.51|) where TZri and 
Erf, dehned by Eqs. (^.49|) , take now the explicit forms 


TZri — - 

^ I+ r] To 


(tZo + r] To TZl 


1 

Ef] = - 

^ l + r]f'o 


{E{no} + vroE{nl}) 


In these equations the quantity tq = Zq/Zq is given by 

ro = [det A/’'tanh —7-fo]^'^^ = J([tanh—7-fi 
2 ^2 


OA 


( 6 . 20 ) 

( 6 . 21 ) 

( 6 . 22 ) 
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where M' stands for the quasi-particle-number operator and T^oa denote the 
eigenvalues of the matrix Tio- The sign of tq is obtained by continuity ac¬ 
cording to (|6.1ti| ). 

The matrix TLq is the solution of the self-consistent Eq. (|5.52|) which, after 
a suitable reorganization, can be written 


Ho {1 — T] tq cotli^ ^Ho) = H . (6.23) 

Using Eq. ( |6.19|) and Eq. ( |5.53|) to dehne H^ from TZq = TZq and H^ from 
TZl, the matrix H is given by 


n 


-r] To 


coth 


H^ 


coth 


(6.24) 


In terms of the matrices TZo and the operator H can be expressed as 






P, 


n = n{no}-'nro j^coth^Tfoj ^{^o} j^coth^Tfoj -2r]rok^ coth , 

(6.25) 


in which the quantity (see Eq. ( b.54| ) and note that k^ = —rjro k^) is 
1 /I Ho " 


k^ = 


l + rjro ( 2 ^ sinh (3Ho 


+ E{nl}-E{no} 


(6.26) 


In this projected variational scheme, Eq. ( p.23D replaces the standard 
HFB self-consistent equation Hq = H{TZo}. The logarithm of the exact N- 
parity projected grand-partition function (Eq. (pTT|)) is approximated by 
InZ^, the stationary value of the functional (|5.46|) , as 

InZ^ = \ti 2 n ln[l -1- e"^^° ] |tr 2 „ [(3HoHri] - pE^j + ln[|(l + T]ro)] . (6.27) 


Equivalently, this quantity can be written as 


\nZ,, = \nZ-T]ro/3k^+\n[^{l + T]ro)] , (6.28) 

where InZ is the logarithm of the HFB partition function, S'{77o} — l3E{TZo}, 
which is associated with the independent-quasi-particle contraction matrix 
Uo (Eq. (^)). 

The entropy Srj of the exact projected state is approximated by the spe¬ 
cialization to the present case of the quantity Sr, given in Eq. (|5.59|) : 

Sr, = -ltr 2 n[H^lnRo + {1 - 1Zr,)ln{l - TZo)] -h ln[i(l-h r^ro)] . (6.29) 
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After some rearrangement, S„ takes the form 


Sr, = 5{7^o} - 


r/ro 


tr 2 n 




2 {l + r]ro) sinhf3Ho 


+ \n[\{l+ T]ro)] 


(6.30) 


where S'jT^o} is the nsnal HFB entropy given by Eq.( |4.6|) . As stressed in 
Sect. 5.4.2, the eqnations (|6.23H6.26|) satished by Hq ensnre that the approxi¬ 
mate self-consistent projected entropy Sr, obeys the standard thermodynamic 
relations (|2.28|) and ( p.29| ) with Fq = —f3~^\nZr,. 

Since onr trial operators F’(m) and A{u) satisfy the validity conditions 
of Eq. (|2.23|) , the density matrix which optimizes the partition fnnction can 
be used to obtain the variational approximation for the expectation value 


of any single-quasi-particle operator Q of the type (|3.12|) . For instance, the 
transposition of Eq. (p.23|) gives : 


{Q) =(1+ QT^r, 


(6.31) 


which generalizes the HFB formula (|4.9|) . Moreover, since these operators Q 
commute with (—1)^, their characteristic function can be found by means of 
the formalism above by changing K into K + ^Q//3, as indicated in Sect. 2.3. 

Altogether, to obtain the variational approximation for the A^-parity 
projected grand-partition function dehned by Eq. O. one hrst deter¬ 
mines the matrix TCq (jointly with the contraction matrix TZr,) from the self- 
consistent equation (|6.23|) which contains (through 7i) the c-numbers ro and 
k^, given by (|6.22|) and (|6.26|) , respectively. One then evaluates Z^ by means 
of Eq. ( 6.28 ). Approximations for the thermodynamic quantities such as 
the entropy (|6.29H6.30|) , or for the expectation values (|6.31|) of single-particle 
observables, in particular of N, involve the averaged contraction matrix TZr, 
dehned in ( |6.20|) ; this matrix is positive since, from Eqs. ( |6.18|) , (|6.19|) and 
(|6.22|) , its diagonal eigenvalues are given by 


FriX — 


l-vYi fanh^TYo/x 


"" e/^Wo. + i B 

1 + T][[ tanh-7do^ 


(6.32) 


According to the general discussion of Sect. 2.3, our variational approxi¬ 
mation for the huctuation of N is given by the thermodynamic relation 


19W 1 , 


(6.33) 
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The derivative dlZr,/dfi is deduced from dTCo /9/r, itself evaluated from Eq. (|6.23|) . 

For large systems, tq is exponentially small while is proportional to the 
volume, as shown by Eqs. (|6.22| ) and (|6.26| ), and these behaviours imply that 
the projection does not modify the HFB results. However, for sufficiently 
small systems and at sufficiently low temperature, the factors entering tq 
(E q. (|6.22|) ) may be close to ±1, except near the Fermi surface ; the terms 
involving tq may thus be significant compared to the dominant terms, al¬ 
though |ro| is always less than 1. 

The results of this Section could have been obtained directly from the 
formalism developed in Sect. 6.1 by using the specific u-dependence of 
and the commutation of Hq with H. The latter property, already established 
in the Remark of Sect.5.4.2 for a more general case, can also be proven 
from Eqs. (|6.14|) and (|6.15|) once the n-independence (|6.18|) of the contraction 
matrix is acknowledged. 


6.3 The Projected BCS model : Generalities 

We now use the formalism of the previous Section by specializing to the 
BCS-type pairing. In this case, the n creation operators a\ can be ar¬ 
ranged (possibly by means of a unitary transformation) in a set of couples 
{(aj,, Op) , 1 < p < 'n./2} such that the Hamiltonian (|3.27|) takes the form 

H ^ ( Cpi^dpdp “ 1 “ dpdpj ^ ) Gpqdj^pdjpdqdq . ( 6 . 34 ) 

p pg 


This may encompass models for electrons in a superconductor (p then denotes 
the momentum and a spin pointing upwards, while p denotes the opposite 
momentum and spin), or for nucleons in a deformed nuclear shell (the mo¬ 
mentum is then replaced by the component of the angular momentum along 
the principal axis of inertia). For the sake of simplicity we have not intro¬ 
duced a magnetic held, which would lift the degeneracy between the states 
p and p. We have also assumed that the diagonal matrix elements Gpp are 
zero (a non-zero value of Gpp entails a small state-dependent shift of the 
single-particle energy e^, which leads to more complicated formulas but does 
not modify the conclusions). Moreover, we consider the usual case where the 
quantities Gpq are real and symmetric in the exchange p q. 

(The reader who wishes to begin with this Section 6.3 without having 
followed the derivations of Sects. 5, 6.1 and 6.2 can do so if he is ready to 
adopt the dehnitions and accept the expressions and equations that we now 
recapitulate. Our variational space was introduced in Sect. 6.2, Eq. (|6.17D , 
where the projection was given in Eq. (|6.2|) and where the trial operator 
Ho (see Eq. (p.44|)) involved a c-number ho (which will play no role in the 
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following Sections) and a 2n x 2n hermitian matrix Tio- The contraction 
matrices TZq and were defined in Eqs. ( |6.18| - ^T^) from the matrix Hq 
which is itself determined self-consistently by Eqs. ( |6.23|) and ( |6.25| ). The 
matrix IZr, and the scalar qnantities tq, and were given in Eqs. (|6.20|) , 
(|6.22|) , (|6.21|) and ( |6.26|) , respectively. The approximate projected partition 
function entropy Sr/ and expectation valnes (Q) were given in Eqs. (|6.28|) , 
dp^) and (lOlD .) 


6.3.1 The Variational Space 

Within the BCS scheme the operator Hq can be pnt into the familiar diagonal 
qnasi-particle form 


Ho = (ho - P Cp) + P ep( 6j6p + bpp ) , (6.35) 

P P 

by means of the canonical transformation 

bp 
bp 


Updp ~\~ Vpdp 


Updp 


Vpdp , 


p Vjpbp Vpbp , 

V Vpbp Vj}b^ 


(6.36) 


p^p 


(Note that two different sign conventions are cnrrently nsed in the literatnre 
for this transformation.) The parameters Up and Vp can be chosen real and 
non-negative, with Up"^ + Vp^ = 1. The variational qnantities are thns e^, Vp 
and ho- 

All the 2n x 2n matrices that enter onr formalism can then be decomposed 
into a nnmber n/2 of 2 x 2 symmetric matrices labelled by p and corresponding 
to the operators dp and a| , pins another set of matrices labelled by p which 
differs from the previons one throngh a change in sign of Vp. In particnlar 
the p-components of the matrix TLq (Eq. (|5.44|) ) are 

Hop = epUp , (6.37) 


where Up is the real symmetric orthogonal matrix 


/ Up Vp 2upVp 

y 2upVp Vp^ — Up^ 


(6.38) 


Accordingly, the p-components of the contraction matrix TZq (Eq. ( |6.18| )) 
read 


H-op — 


Ti 


Trfl'^](;5) 


dpdp 

dpdp 


Qjp Qjp 


= l{l-Uptp) 


(6.39) 
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where we have introduced for shorthand the notation 


tp = tanh|/?ep . 

Likewise, the matrix TZ\ (Eg. ( |t).l9| )) has the p-components 


(6.40) 


(6.41) 


The energy E{71q], associated with the operator K = H — fiN and 
evaluated for ( |6.39|) , is now given by 

^ ) i,^p /^) [1 i.^p ^p )^p] ^ ) ^pq^p^p^q'^q^p^q • (6.42) 

p pq 

The corresponding matrix TipjT^o} has the p-components 


'Hp{TZo} — 


^p-^^ ^P 


with 


/\p — ^ ) ^pq ^q^q ^q 


(6.43) 


(6.44) 


By substituting ^ to tg in Eqs. (|6.42|) and ( |6.43|) we obtain the energy 
E{7ll} and the matrix In particular, is replaced by 


Ap — ^(0 Gpg UqVq 


-1 

q 


(6.45) 


The c-numbers pq (Eg. |6.22|) ) and k} (Eg. (|6.26|) ) are respectively equal to 

ro = n ^p^ ’ (6-46) 


and 


= 


1 + pro 


(^p ^p) { ^p i^p k) i'^p '^p ) UpVp (Ap + Ap) I. 

(6.47) 


6.3.2 The Variational Equations 


p and Vp must be determined from the self-consistent 


The trial parameters e 
matrix equations (|6.23|) and (|6.25|) . By means of the canonical transformation 
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( |6.36| ), let us write this equation in the basis where Tio is diagonal. From the 
diagonal element of Eq. (|6.23|) we thus obtain: 


ep{l-rirotp ‘^) +27]= 

(Mp^-n/)(ep-/i)(l-?7rotp“^) + 2npnp(A°-?7roAjtp“^) , ‘ ^ 

and from the off-diagonal element : 

0 = 2npnp(ep -/r)(l-h r/rofp"^) - (wp^ - V)(Aj-F? 7 ro Ajfp-^) . (6.49) 


Together with the preceding definitions of A®, A^, fp, tq and k^, Eqs. (|6.48|) 
and ( 6.49|) are the self-consistent equations which determine the variational 
parameters Cp, Up and Vp. 

Let us further work out these equations. By introducing the weighted 
averaged pairing gap 


^ Ag + ?/roA^tp ^ 
1 + Tirotp~'^ 


(6.50) 


and making use of the relation {up^ — Vp^Y + {2upVpY 
Eq. (|6.49D into 

2upVp Up^ — Vp^ 1 

A77 ri Cp ]1 Cp 


VP 

where Cp has the familiar BCS form 


Cp — 


\lYp ~ dY + 


1 , we can simplify 

(6.51) 

(6.52) 


Except for the definition of A^p, Eq. (|6.51| ) is the same as the usual BCS 
equation for the coefficients Wp, Vp of the canonical transformation (|6.36|) . 
As regards Eq. (|6.48|) , it appears as an extension of the BCS relation Cp = 
(up^ — Vp^){ep — ]i) + 2upVpAp ; actually, we can rewrite Eq. 


as 


Cp = (Mp2 - VpY (Cp - /i) + 2UpVp[App - -2 -‘2 (- A°)] 


2r] Pq k^ 


f ^ _ V^f 


(6.53) 


tn-prot^^ 


'^p ■/ ■ U “p 

and hence, by using Eq. (|6.51| ) to eliminate Up and Vp, as 

2prok^ 2pro App(Aj - A 


Cp Cp 


tp-pro t. 


-1 

p 


Lt) > nt'- 


2+ -2 


(6.54) 


0''P 


These quantities Cp are the energies which appear in the diagonal form of 
the density operator. Eqs. ( |6.52|) and ( |6.54|) show that Cp has no longer 
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play here the 
and 


the BCS form. Indeed, several quantities, A^, and A^p, 
role of pairing gaps. We recover the standard BCS relation between Cp 
App when the system becomes sufficiently large, in which case tq and tq 
become small. We can eliminate Up and Vp from the dehnitions (|6.44|) , (|6.45|) 
and (16.501) of A^, A^ and A^p. These quantities thus satisfy the equations 


AJ-AS=2 


1 ^ Apg 

<? 


(tq ^ - tq) 


(6.55) 


^vp “ 2 ^ 


A, 




tq l + T] To tp ^tq ^ 


-2 


(6.56) 


q eq l + r]rotp 

Eq. (|6.56|) exhibits the projection correction to the so-called BCS gap equa¬ 


tion, while Eq. ( 6.55 ) obviously has no BCS counterpart. Finally, by using 
(|6.53|) and the identity Y.p '^p'>^p{‘^p'tp~ ^ptp^) = 0, we can write the expres¬ 
sion (p.47]) for in a more explicit form : 


1 + 


2r] ro 


E 


tp-^ -1 


l + r]rop l-r]rot 


-2 


E 


tp tp 1 -\- rjPQtp ^ App 


(a; - A5) 


2{l + r]ro) 1 -r^rotj 


-2 


(6.57) 


6.3.3 Comments on the Equations 


Altogether then, we have thus to solve the self-consistent equations (|6.54|) 
and ( |6.56| ) for the two sets Cp and A^p. The quantity Cp is indirectly involved 
through tp = tanh(/dep/2). The two c-numbers tq and are expressed by 
(|6.46|) and (|6.57|) , respectively. The difference A^ — A^, which enters (|6.54|) 
and (|6.57|) , can be eliminated by means of Eq. (|6.55|) . The coefficients Up 
and Vp are given explicitly by ( |6.51| ). 

In Eq. ( |6.56|) the last fraction, which takes care of the projection effects, 
can also be written as 1 -|- (tq~‘^ — 1) f {rjrQtp~‘^) with /(x) = x/(l -|- x). The 
positive term {tq~‘^ — 1 ) differs signihcantly from zero only when (3eq < 1 . 
For an even-number projection {rj = 1), the function / is positive and varies 
between 0 and 1/2 while for an odd-number projection (—1 < r]ro < 0) / is 
negative and varies from — cx) to zero. The following conclusions can therefore 
be drawn : 


in the sums of Eqs. (|6.56| - |6)57D , the modihcations due to projection af¬ 
fect predominantly the terms corresponding to the single-particle states 
near the Fermi surface. 
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• the even-number projected gap is larger than the BCS one while the 
odd-number gap is smaller. 

• larger differences with BCS are expected for the projection on an odd 
number of particles. 

Once Eqs. (|6.54 ) and (|6.56 ) are solved, the thermodynamics is governed 
by the equations (|6.30| ), (|6.20D and ( |t).23| ) which give 5^, and re¬ 
spectively. To compute the approximation Srj given by (|6.3CI|) to the exact 
number-parity projected entropy 5^, we take into account the twofold degen¬ 
eracy of the eigenvalues (e^^®p -|- 1)“^ of the contraction matrix TZq. We find 
that Sr, differs from the independent quasi-particle entropy S'bcs, 


-Sbcs = S {7^o} = cosh - Pcp tanh , 

p / / 


according to 


Sr, — SbCS — 


r]ro 


1 , E/^ep(tp^-tp)+ ln[|(l + 77ro)] 

i- + p 


(6.58) 


(6.59) 


It results immediately from ( |6.31|) , ( p.39| ) and ( |6.51| ) that the expectation 
value of the particle number is given in terms of the chemical potential by 


{N)=T. 


p I 


6p — tp Tj T^t 


-n 


1+ r]rQ 


(6.60) 


This expression is also the derivative with respect to /i of Fq = —\n.Zr,/j3 = 
Er, — Sr,/13, which is our approximation to the logarithm of the exact N- 
parity projected grand-potential —InZ^/ /3. The expectation value {H — fiN) 
is obtained from Eqs. (|6.21|) and (|6.42|) as 


{H-^N)^Er,= Y.i^p-^^) 


[U, 


2 2^^P + ^^0^p 


-n 


v/) 


I + 7] To 
-1 .-1 


^ ) GpqUpVp UpVq 


tptq + 7] Tq tp tg 


(6.61) 


pq l + vro 

while, using Eqs. ( |6.5(J| ), (|6.51|) , ( |6.60[ ) and ( |6.61| ), the energy {H) is given by 

^ ^ ^ I -t- n l-r^l 

(6.62) 




p i 


Cp(cp - /^) + 2\p ^p + Vrot 


-1 

p 


1 +7] Tq 


Here again, the deviation from the BCS energy is exhibited by the two terms 
that contain tq. 
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In accordance with the general properties of our variational scheme, the 
thermodynamic relations ( p.27|) and (|2.28|) are satished by the approximate 
entropy The expectation values of the single-particle observables are 
given by ( |6.31|) . The specihc heat, as expected, is equal both to —/9d5^/d/d 
and to the derivative of {H) with respect to the temperature. In Eq. (|6.62| ), 
the temperature appears directly in tp = tanh(/9ep/2), and indirectly through 
Cp, App, ro and also through fi since (N) should be kept fixed in the derivative. 
Finally, the general arguments at the end of Sect.6.2 show that the charac¬ 
teristic function for any single-quasi-particle operator commuting with K, 
such as N, is variationally obtained by adding a source term to the operator 
K. We hnd therefore that the fluctuation satishes 


AN^ = 


1 d{N) 
f3 d/i 


(6.63) 


where in Eq. (|6.6CI|) the derivative should be taken both explicitly and through 
Cp, App and tq. 

A somewhat similar extension of the BCS theory, also involving a pro¬ 
jection on the spaces with a given parity of the particle number, has been 
worked out in Ref. ||16||. The results agree with ours for the partition func¬ 
tion Zp. The comparison between our variational equations (|6.48|-|6~i9|) and 


those of Ref. [16 


why the entropy found in Ref. 


appears more difficult. 

TBI 


In particular it is not clear to us 
does not satisfy, in contrast to ours, the 


thermodynamic relations ( 2.27|) and ( 2.28 ). On the other hand, the general 
discussion of Sect. 2 as well as Eq. (|6.63|) show that fluctuations and corre¬ 
lations calculated with our method must definitely be different; for instance 
we find here that the fluctuations of the particle number tend to zero with 
the temperature, as expected from general grounds. 


6.4 The Projected BCS Model : Limiting Cases 

6.4.1 Low Temperatures 

At zero temperature the BCS (and HFB) solutions of the ground and excited 
states are exact eigenstates of the number-parity operator and their eigen¬ 
values ±1, depend on the number of quasi-particles added to the even BCS 
vacuum. 

Let us first recall the low-temperature limits of some quantities and equa¬ 
tions of interest. For large values of /?, the usual BCS gap equation becomes 

ABCSp=iEG'„^(l-2e-'3'») . (6.64) 
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In Eq. ( |6.64|) the terms proportional to exp {—f3eq) lead to a decrease of the 
gaps versus temperature T = 1/(3. In the same limit, the BCS entropy (|6.58|) 
is approximated by 

^Bcs = 2 5:/3epe-^''p . (6.65) 


In the physics of mesoscopic superconductors, one often considers a sim- 
plihed situation in which the twofold degenerate unperturbed energies Cp are 
sufficiently dense near the Fermi surface to allow the replacement of a sum 
over p by an integral over e = ep — p with a weight wp equal to the level 
density. One also assumes that the pairing matrix elements are constant and 
equal to G over an interval A around the Fermi energy p and that they vanish 
outside. The width of this interval, taken of the order of the Debye energy, is 
assumed to be much larger than the pairing gap. It must be hnite, however, 
to ensure the convergence of the integrals. The ordinary BCS gap equation 
at zero temperature 

2 M/2 wp de 

G ~ i-A/2 + A2 

dehnes a p-independent quantity A that we shall take below as a reference 
energy. In this held of physics, the “low-temperature” regime actually refers 
to an intermediate range of temperatures such that the conditions 



i -C ^ -C ^ < wf (6.67) 

are valid. Under these conditions Eqs. (|6.64D and ( |6.65[) give, after elimina¬ 
tion of the cut-off A by means of Eq. (|6.66|) , analytical expressions for the 
“low temperature” dependence of the BCS pairing gap : 


Abcs(/5) = A 



( 6 . 68 ) 


and of the BCS entropy : 

Ages = 2 wfA Y^27r/3Ae“^^ . (6.69) 

In the A-parity projected extension, the quantity tq (Eq. ( |6.46|) ) plays 
a central role. In the “low-temperature” limit, its value tends in general 
towards 1 and one has 


Inro = —4 ^ e 

p 


(6.70) 
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up to terms equal to, or smaller than exp(—3/9A). (However, in the case 
wpA ~ 1 that we shall numerically investigate in Section 6.5.4, the condi¬ 
tion (|6.67D cannot be satished and we shall see that tq vanishes with the 
temperature.) 

We shall check below (Sect. 6.5.2) that, in the limit wpA 3> 1, the pairing 
gap Arip obtained with the iV-parity projection does not depend much on p 
and that it is almost identical with Abcs(/ 5)- This will be illustrated by the 
numerical results plotted i n Fig. 2. Moreov er, we shall see that the difference 

ep-/i)2 +A2^ ~ + A2 is small. Under the 


between Cp and e(e) = 


VP 


conditions (|6.67|), the approximation 






8A3' 


/27rA 

~T 


He) 1 + 




(6.71) 


holds. Hence the expression (|6.70|) becomes 


In To = —4 wp 


1271A 


-/9A 


/3 


1 + 


8pA 


(6.72) 


We can then evaluate the odd-even grand potential difference for given 
values of T and /i, 

6FG = ~{lnZ_-lnZ+) , (6.73) 


with good accuracy in terms of the BCS solution. From Eqs. (|6.28|) and 
(|6.72|) , neglecting the contribution ro{kf -|- k\), we obtain for T in the range 

(CT) 


^ ^1 ^ (87r(M;pA)2 ^ 

A (3A 1 + ro 2A V ^ ^ A 


(6.74) 


Ref. 1^ reports on a direct measurement of the difference SF between the 
free energies of neighbouring odd and even systems for small superconducting 
metallic islands and shows that it agrees remarkably well with the approxi¬ 
mation ( 6.74 ). We shall return to this point, and to the difference between 
5Fq and 5F in Sect. 6.5.3. 

To proceed further, and in particular to deal with very low temperatures 


(3 ^ wf, for which (|6.67|) no longer holds, we consider separately the pro¬ 
jections on an even or odd number of particles. Let us begin with the even 


case. 
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Even Particle Number 

Using (|6.70|) , the expansion of Eq. (|6.56|) yields: 


9 {<l¥^p) ‘ 


M _ 2e — 4 e | 

V {r^P,<l) / 


(6.75) 

In the low-temperature limit (3/wf 3> 1, Eqs. ( |6.64|) and (|6.75|) are identical. 
In Eq. (|6.75|) , the terms which determine the low-temperature behaviour of 
the gaps are negative. However, they are now of the order of, or smaller than, 
exp (—2/3A) while they were of the order of exp {—f3A) for BCS. This reflects 
the fact that any excitation with hxed iV-parity requires the creation of at 
least two quasi-particles and therefore an energy larger than 2A (instead 
of A for BCS). Therefore the decrease of the pairing gaps with increasing 
temperature will be slower than for the BCS case. The low temperature 
behaviour of the entropy (|6.59|) is now given by 


= 2 ^ f3ep e“2^ep ^2 ^ f5{ep + ej + ^ 9 ) 

P pq (p^q) 


(6.76) 


It is also governed by terms smaller than exp (—2/dA). Under the conditions 
(|6.67|) , the gap equation (|6.75|) provides an analytical expression for the low- 
temperature dependence of the positive-parity-number gap A+(/9) : 


A+(/3) = A 


l_4vr^e-2/^^ 


(6.77) 


The extrapolation to zero temperature of A+ is therefore equal to the BCS 


gap A. Under the same conditions, one can can rewrite the expression (|6.76|) 
as 

= Svr (tcpA)^ e~2/^A ^ (6.78) 

The BCS entropy (|6.69|) is therefore larger than Indeed, the disorder 
decreases when odd-particle components are removed. 

Terms proportional to exp (—2/3A) occur in other thermodynamic quan¬ 
tities. For instance, for the average particle number, the limit of large (3 

is 


(A) = E 


/i 


1 - 2e 


-2/?e, 


4 E e 

9 


(3(ep T Cq) 


(6.79) 

The analysis of Eqs. ( |6.54D , (|6.55|) and (|6.57|) shows that, in the large [3 
limit, the energy Cp of a state in the vicinity of the Fermi surface is smaller 
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than its BCS counterpart. The correction depends on the two degenerate 
quasi-particle states with lowest energy. Hereafter we shall denote Aq and Cq 
the corresponding gap and energy. Assuming moreover that all the matrix 
elements Gpq (p 7 ^ q) have the same value G, one hnds 



(6.80) 


Actually, in the numerical analysis of Sect. 6.5, we shall hnd that and Cp 
differ by only few percent. 

Odd Particle Number 

Let us now study the case of an odd system. It is useful to compare our 
results with those obtained from two different types of BCS approximation. 
The first one (hereafter labeled BCS) corresponds to the standard equation 
with the usual constraint on the odd average number of particles. The second 
one (labeled BCSl) corresponds to a situation in which there is exactly one 
particle on the twofold-degenerate individual level whose energy e is closest 
to the Fermi energy qi. If we label this level with the index p = 0, the 
BCSl solution is obtained by effecting in all sums the substitution (uq, Vq) —>• 
(uo, —Uq) for only one of the two quasi-particles states 0 or 0 associated 
with cq. Under this transformation the product uqVo changes sign and the 
contributions of these two quasi-particles states to the gap cancel exactly. 
Thus for p 7 ^ 0, the BCSl gap equation at arbitrary temperatures is 



(6.81) 


«(«#!)) 


The only difference from the usual BCS equation is the absence of the g = 0 
term in the sum. However, through self-consistency, this affects the value 
of the gaps. (A constraint on the odd value of the average total number 
of particles, which includes the single particle occupying the level 0 , must 
of course be enforced.) In the BCSl solution, the particle which occupies 
the state p = 0 (or p = 0) forbids the formation of a pair (0, O). Within the 
single-particle space from which the levels 0 and O are omitted, the variational 
pairing equations must then be solved for the (even) system with one fewer 
particle. We shall see (Fig. 1) that this results in a sizeable difference between 
BCS and BCSl. 

Returning to the odd-number projection, we have to treat separately the 
gaps A_p for the quasi-particles p 7 ^ 0 and the gap A_o for p = 0. Let us 
first consider the case p 7 ^ 0 . At low temperature both the numerator and 
the denominator in the fraction in the right hand side of Eq. (|6.56|) vanish 
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as exp(—/?eo). After dividing out by this common factor, we obtain the 
well-conditioned equation: 


^-p-n \^PQ 

q{q^0) I 


A, 


+ a 


A, 


pO" 




A_ 


Co 


pq 


l] - eo) 


(6.82) 

The zero-temperature limit of Eq. (|6.82|) differs from that of BCS, but is 
identical with the limit of ( |6.81| ). This implies that, as > oo, Eq. ( |6.56| ) 
automatically singles out the level with the lowest quasi-particle energy and 
generates the odd-particle-number ground state given by BCSl. The low- 
temperature correction to the gap A_p is positive and depends on terms 
which are proportional to an exponential of the difference between a quasi¬ 
particle energy Cq and cq. Therefore, starting from the value provided by 
the BCSl equation for > cx), one expects the gap to hrst increase with 
temperature as exp (—/9/2 i<;f^A). 

The case p = 0 is slightly more complicated. One hnds that the temper¬ 
ature corrections involve no longer cq but the second smallest quasi-particle 
energy ei, the corresponding gap A_i and degeneracy gi : 



A-n- ■ „ _ 

2 ei ^ 

,-f3{eq - ei)| 
(6.83) 

Despite their formally different expressions, we shall find in the numerical 
application performed in Sect. 6.5 that the values of A_o and A_p are equal 
within a few percent. In nuclear physics, one often encounters situations 
where 1/A, /3 and wp are of the same order of magnitude. Then, the very 
low-temperature entropy (|6.59 ) behaves as 


^_ = ln2 + ^ (3 (cp - eo) e“^(^P “ . (6.84) 

p (p^o) 


When (3 ^ oo, the entropy is equal to ln2, a value which reflects the 
twofold degeneracy of the ground state, and to which we shall return in 
Sect. 6.5.4. 


Under the conditions (|6.67|) , using Eqs. (|6.71|) and (|6.72|) , one can refor¬ 
mulate Eq. (|6.56|) as an equation for the (p-independent) odd-number gap 
A_(/5) : 


2 _ f^/‘^ wpde _ 1 ^ 1 

G 2-A/2 ^e^ + A_{f3y ^ 2/3A_(/3)2 
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The limit of this equation, in the regime (|6.67|), yields 


A_(/?) = A 


1 - 


+ 


2wfA d/dtapA^ 


( 6 . 86 ) 


At zero temperature the odd-A gap A_ is therefore smaller than the BCS 
gap by the quantity 1/2wf. This point has already been mentioned in con¬ 
nection with the discussion of the low-temperature gap equation (|6.82|) and 
its relation with the BCSl Eq. ( |6.81|) . In the temperature range defined by 
the conditions (|6.67|) , the gap grows linearly with temperature. We have 


already noted that it grows as exp{—P/2wF‘^A) when [3 is larger than wp, 
or of the same order. The odd-A projected thermodynamic quantities have 
a larger temperature variation than the corresponding BCS ones. This re¬ 
sults from the fact that excitation energies are then of the order tcp”^ and 
therefore much smaller than A . 

For the entropy (p.59|) and for low temperatures satisfying the conditions 


(|6.67|) , we obtain : 


^_ = ln2 + - + -In 
2 2 


^TTtCp^A' 


(6.87) 


This is simply the expression of the Sakur-Tetrode entropy m for one clas¬ 
sical particle of mass A with an internal twofold degeneracy, enclosed in a 
one-dimensional box of length 27rhtcp. Such a result is natural since the 
elementary excitations concern only a single unpaired particle, with energy 
Cp ~ A -|- e^/2A; if we identify e = Cp — with the momentum of a particle 
in a box, the identification of the level densities leads to the size 2nh wy for 
this box. While elementary excitations of even-A systems require energies 
of at least 2A, the effective spectrum for odd-A systems is nearly contin¬ 
uous, with the form m^/2r(;|A where m is an integer. This occurence of 
gapless elementary excitations is consistent with recent experiments on the 
conductance of odd-A ultrasmall Aluminium grains[^]. The analogy with 


the Sakur-Tetrode entropy shows also that in the variational projected solu¬ 
tion, contrary to the BCSl one, the unpaired particle does not remain on the 
single-particle level of energy cq = p. The projected solution is a coherent su¬ 
perposition of conhgurations in which this particle explores all single-particle 
levels consistent with the temperature. The fact that we recover a formula 
of classical thermodynamics is another indication that, strictly speaking, the 
regime associated with the conditions (|6.67|) is not a low temperature regime. 
Actually, ( |6.67|) implies that the temperature is sufficiently large so that 
wy^A/(3 S> 1, which ensures the positivity of the entropy as given by 
Eq. (^:h7D . 
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For the average number of particles^ the low temperature limit of Eq. ( |6.6(]|) 
yields 


(Af) = i+ y 

p ip¥=o) 




1 -2e' 


-(3{ep - eo) 


( 6 . 88 ) 


In contrast to the corresponding BCS formula, the temperature correction 
depends also on the difference Cp — Cq. 

Using Eqs. ( |6.54D , (|6.55| ) and ( |6.57[ ) one Ends, in the limit jS/w^ 3> 1, 


that Op is larger than Cp. Taking again all the matrix elements Gpq {p ^ q) 
equal to the same value G, one has 


_ G 

Cp — Cp + — 






ei 


(6.89) 


To derive this formula, we have taken = 2 as in the numerical examples 
presented in Sect. 6.5. In these applications, we shall hnd that the magnitude 
of Cp — hp is small (a few percent). 

The inequalities S- > S'bcs > implied by Eqs. ( |6.69| ), ( |6.78|) , (|6.84|) 


and ( |6.87| ) are easily understood. When the system is constrained to have an 
even particle number, we have seen that its lowest excited states are obtained 
from the ground state by the creation of two quasi-particles, which requires 
at least an energy 2A. For an odd system, the excitation of the (gapless) 
unpaired last particle is much easier. At non-zero temperature the BCS state 
is a weighted sum over even and odd components. The concavity of entropy 
implies that S'bcs is larger than the smallest of the quantities and 
that is S+. The inequality S'bcs < is made possible by the very small 
weight of the odd components entering the BCS state at low temperatures. 

The smallness of this weight also explains why, as already noted in Sect. 6.3.3, 
the odd-A projection affects the BCS results more than the even one. Indeed, 
we have already met with several qualitative differences between the prop¬ 
erties of odd-systems and their BCS description, while BCS explains fairly 
well the even systems. This point will also be exhibited by the numerical 
applications of Sect. 6.5, especially on Figs. 1 and 6 which show the gaps at 
the Fermi surface as a function of temperature. 


6.4.2 High Temperatures 

High temperatures relevant for superconductivity are of the order of the criti¬ 
cal temperature, that is Tc = 0.57A. General conclusions on the behaviour of 
the equations can only be obtained when the density of single-particle states 
at the Fermi surface is significantly larger than one, which entails wpA S> 1. 
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In this limit, the quantity llr^p.g tends towards 0 for any value of p and 
q. Then, Eq. ( |6.56|) is equivalent to : 


^VP=\ Y. Gpq^tq[i+p{l-tq^) H ' (6-90) 

^ gi<lT^P) I rytp,q ) 

This differs from the BCS equation by the factor within the braces, which 
deviates from 1 by a vanishingly small correction. The projected gap values 
converge (from above when r; = 1 or from below when rj = —1) towards 
the BCS values as one reaches the critical temperature. In particular, the 
odd solution converges towards the standard BCS solution and not towards 
BCSl. 

In Sect. 6.5.4, we present a numerical case which does not pertain to the 
limit 3> 1 and yields qualitatively different results. 


6.4.3 Large Systems 


In large systems, at non-zero temperatures, the product ro decreases as the 
inverse of an exponential of the volume of the system. Then, an iterative 
solution of Eqs. (|6.54|) , (|6.56|) , (|6.62|) and (|6.57|) is attainable by starting 
from the BCS solution. Since in several places pq appears multiplied by 
coth^(/3ep/2), the effect of the projection is surmised to be especially sizeable 
near the Fermi surface. 

In this large system limit, the approximate projected entropy (|6.59|) , for 
both values of rj, differs from the BCS entropy by a constant: 


Sp = Sbcs — ln2 


(6.91) 


This result, shows that in large systems (for a given temperature) the num¬ 
bers of states with odd and even A-parity are equal. 

The quantity increases linearly with the volume, so that Cp — Cp is ex¬ 
pected to be dominated by the second term in the r.h.s. of Eq. ( |6.54|) . Then, 
at the lowest non-trivial order, the quasi-particle energy can be expressed as 


~ e„ + 




qq' 


G, 


qq'- 


Apg (t tq) A.pq> {t , tqi) 


€q' 


(6.92) 


The self-consistency and the constraint on the particle-number moreover shift 
fi and App away from their BCS values. 





















6.5 The Projected BCS Model : A Numerical Illustra¬ 
tion 

To evaluate the corrections to BCS introduced by the number-parity projec¬ 
tion, we have performed three schematic calculations that exemplify situa¬ 
tions encountered in the physics of (i) superdeformed nuclear states around 
mass A = 190, (ii) superconducting mesoscopic metallic islands and (hi) 
ultrasmall metallic grains. 


6.5.1 Superdeformed Heavy Nuclei 


We hrst consider an application to a held of nuclear structure physics which, 
recently, has been the subject of an intense experimental exploration. Su¬ 
perdeformed nuclear states are generally populated via fusion reactions be¬ 
tween two heavy ions. Once the colliding nuclei have fused and a few neutrons 
have evaporated, the compound nucleus decays with a small probability (a 
few percents) into superdeformed metastable conhgurations via the emission 
of statistical y-rays. In heavy nuclei, the angular momentum dependence 
of the moment of inertia observed in superdeformed bands gives strong evi¬ 
dence of the presence of pairing correlations. Most of the presently available 
data concern the properties of superdeformed ground states or of the lowest 
excitations. However, it is expected that the next generation of y-arrays, 
in Europe and in the US, will soon allow investigations of the statistical 
behaviour of the spectrum and of the properties of heated nuclei. 

In view of the general scope of the present work, rather than studying a 
realistic nuclear conhguration (see for instance Ref. [Q) that would require 
the solution of the full extended HFB equations written in Sect. 6.2, we 
prefer to investigate a simple model in which the nuclear context outlined 
above is used as a guidance for selecting the relevant parameter values. Thus 
we consider twofold degenerate single-particle levels of energies e^, regularly 
spaced with a density wp =2MeV“^. Our active pairing space includes 21 
levels and extends 5 MeV above and 5 MeV below the Fermi surface, so 
that the cut-off A equals 10 MeV. The matrix elements of the interaction are 
taken to have the form Gpg = G (1 — 6pq) with G = 0.23MeV. For this value 
of G, the BCS pairing gap A at zero temperature is IMeV, consistent with 
the data for the even heavy nuclei. In fact, it is partly because A is close to 
IMeV in medium and heavy nuclei that the MeV turns out to be a natural 
energy unit in nuclear dynamics despite the fact that binding energies are 
several orders of magnitude larger. 

For these values of the parameters wp, A and G (or equivalently A), 
we have solved the A-parity projected coupled equations ( |6.54j6.57| ) with 
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r] = ±1 as function of the temperature. This provides us with the projected 
pairing gap A^p defined by Eq. (|6.5CI|) . Although our interaction is very sim¬ 
ple, it leads to state-dependent gaps for the BCS as well as for the projected 
cases. The numerics, however, shows that this state dependence is weak: we 
hnd that at all temperatures the variation of the gap A^p versus p never ex¬ 
ceeds 9%. We thus content ourselves with plots of Ap, that is, of the value of 
the gap associated with the quasi-particle having the lowest energy in either 
the even or odd system. 

The quantity A which serves as a reference energy in all the figures of 
the present Section, both for the temperature in abscissa and for the gap Ap 
in ordinates, is likewise the zero-temperature BCS gap associated with the 
lowest energy quasi-particle. 

We hrst consider (upper part of Fig. 1) the even-number ca.se with (N) = 
20. (Here, (N) should not to be confused with the numbers of neutrons and 
protons which, for a heavy nucleus, are of the order of 120 and 80, respec¬ 
tively. (N) corresponds to the particles within the active pairing space; it 
does not take into account the nucleons of the fully occupied levels whose 
energies are smaller than p — A/2.) We have plotted the temperature de¬ 
pendences of the BCS gap and of the projected gap Ap. As expected from 
the discussion of Eq. (|6.56|) in Sect. 6.3.3 the even-A projected gap is always 
larger than the BCS one. The two curves meet at T = 0 and also at the 
critical temperature Tc = 0.57A, consistent with the analysis of Sect. 6.4.2. 

The lower part of Fig. 1 presents results for the odd-number case with 
(A) =21. It compares the temperature-dependent gaps obtained with our 
odd-A projected method, with BCS and with BCSl. As explained in Sect. 6.4.1, 
in the BCSl option the level p = 0 (with cq = /i = 0) is always occupied 
by exactly one nucleon and therefore not available for pairing. We have 
therefore chosen to show in all three cases the gap Ap associated with the 
lowest-energy quasi-particle whose index is distinct from 0. The odd-A pro¬ 
jected curve always lies between the BCS and BCSl ones. In agreement with 
the discussion of Sect. 6.4.2, the projected curve joins BCSl at T = 0 and 
BCS at the critical temperature. The failure of BCS at low temperature is 
associated with the small weight of the odd-A components in the BCS state. 
The BCSl approximation, by enforcing the odd parity at zero temperature, 
approaches at low temperature the projected state better than does BCS. 
On the other hand, when the temperature increases, the even and odd com¬ 
ponents tend to have equal weights, while the differences between even and 
odd systems fade out. In that limit the BCSl approximation becomes inad¬ 
equate since it prevents the formation of a pair in the level eo. This inhibits 
a correct description of pairing correlations which, at high temperature, take 
place coherently on all levels in the energy range of the gap for the odd as 
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Figure 1: Comparison of BCS and number-parity projected gaps at the 
Fermi surface versus temperature for tcpA = 2, tcpA = 20 and {N) = 20 or 
21. These parameter values are illustrative of the physics of heavy nuclei. 
The upper (respectively lower) part of the hgure corresponds to an even 
(respectively odd) average number of particles. With respect to BCS, the 
pairing gap is enhanced by even-iV projection and (more strongly) reduced 
by odd-iV projection, except near T^. For the BCSl curve, see Sect. 6.4.1 
and 6.5.1. The odd-iV projected gap increases with T up to 0.4A. 
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well as for the even systems. 

When the temperature grows, the odd-A^ projected gap begins to increase 
as expected from (|6.82[ - |OBD , then reaches a maximum, decreases and con¬ 
verges asymptotically towards the BCS curve near the critical temperature. 
At low temperatures, there is a significant difference between the usual BCS 
gap and Ap. In nuclear physics, this difference should be taken into account 
when mass differences between odd and even neighbour nuclei are used to 
extract pairing gap values. 

In the numerical example of Fig. 1, the condition wpA S> 1 is not satisfied. 
It is therefore not surprising that the results do not comply with Eqs. (|6.77|) 
and (|6.86|) derived under the conditions (|6.67|) . 

When gaps and temperatures are measured in A units, the BCS curve is 
universal up to small differences resulting from the details of the spectrum. 
For instance, the critical temperature always comes out equal (with an ac¬ 
curacy better than 1%) to the value 0.57A which is derived for a continuous 
spectrum with constant level density. In contrast, as anticipated in Sect. 6.4, 
the number-parity projected curves show more diversity. In our model, their 
characteristics are mostly determined by the product wpA and to a lesser 
degree by the product wpA and the ratio A/A. The quantity wpA, “the 
number of Cooper pairs”, is the number of twofold degenerate single-particle 
levels in an energy interval equal to the gap, while wpA is the number of 
levels in the active pairing space. In the nuclear physics case of Fig. 1 we 
have taken wpA = 2 and A/A = 10, the number wpA of levels being 21. 
The data corresponding to condensed matter physics are expected to be quite 
different. 


6.5.2 Mesoscopic Metallic Islands 

For superconducting materials, A is of the order of the Debye energy while A 
is of the order of the critical temperature, so that A/A is larger than above. 
On the other hand, depending on the size of the metallic islands or grains, 
such as those considered in recent experiments]^, |], 3, the product wpA 
can vary between 0.5 and several thousands. 

Let us first consider the case wpA S> 1, which corresponds to most of 
the experimental studies performed on superconducting metallic islands. In 
this limit the effects on the pairing gap Ap of the number-parity projection 
are weak. In order to be able to visualize them, we have considered in 
Fig. 2 not too large a value of wpA, taking tcpA = 10 with A/A = 20 and 
hence wpG = 0.34. Thus, the active space comprises 201 twofold-degenerate 
levels. The average particle numbers of the even and odd systems are 200 and 
201, respectively. The BCS and even-A projected curves for Ap are almost 
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indistinguishable. In agreement with the discussion of the projected gap 
equation in Sect. 6.3.3, the odd-A^ projected gap differs from the BCS value 
more than the even-A^. At low temperature, it increases and is represented 
with good accuracy by the linear T-dependence of Eq. (|6.86|) . 
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Figure 2: Comparison of BCS and number-parity projected gaps at the Fermi 
surface versus temperature for wpA = 10, wpA = 200 and {N) = 200 (upper 
part) or 201 (lower part). Such parameter values are illustrative of super¬ 
conducting mesoscopic islands. The effects are qualitatively the same as in 
Fig. 1 but weaker. Even-odd effects disappear above Tx = 0.3A. 


In Fig. 3 we have plotted the entropy of the odd system calculated in four 
different ways : with or without (i.e. G = 0) pairing interaction and with 
or without projection on number parity. The odd-A projected entropy S- 
tends to ln2 when T —0; above the critical temperature it joins the non 
interacting odd-N projected entropy. In contrast, the BCS entropy vanishes 
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at small temperature; above the critical temperature, it joins the grand- 
canonical entropy of the non-interacting system. For T larger than 0.35A 
the odd-iV projection lowers the BCS entropy, according to Sbcs — = ln2. 

In contrast, at small T, one has 3> Sbcs and the projected entropy 
approaches the Sakur-Tetrode regime (Eq. (|6.87 )). Although we are working 
with a discrete spectrnm, for T/A > 0.1, the grand-canonical nnprojected 
entropy of the non-interacting system is almost undistinguishable from the 
expression 

27r^ 

S= - wfT , (6.93) 


obtained by assnming a continuous spectrum with constant level density. 
The only effect of the odd-A projection is again a lowering of the entropy by 

ln2. 

In Fig. 4 we show the specific heat Cy = d{H)/dT = TdS'/dT for the 
odd system. With or withont projection, above the critical temperature 
Tc = 0.57A, the specific heat of the interacting system follows the linear 
form (|6.93|) . Both the critical temperature and the magnitude of the dis¬ 
continuity of Cy are the same for BCS and odd-A projected curves. The 
low temperature behaviour, however, is affected by the projection. The en¬ 
larged scale inset in Fig. 4 shows that the BCS specihc heat drops to zero 
at low temperature, while the odd-A projected curve approaches the value 
I expected for a single classical particle moving in a one-dimensional box 
(see Sect. 6.4.1). As mentioned in Sect. 6.3.3, the odd-A specihc heat Cy 
takes this value | only over the intermediate range of temperatnres where the 
quasi-particle spectrum Cp is well approximated by a quadratic function of 
Cp — pL. In Sect. 6.5.4 we shall investigate in more details its zero-temperatnre 
limit. In the absence of pairing interaction Cy does not depend on whether 
odd-A projection is or is not effected, and it is again given by the linear 
function ( |6.93|) . 


6.5.3 The Free-Energy Difference 

We now consider a quantity which is experimentally accessible: the difference 
5F between the free energies of the odd and even systems [pi|]. Fig. 5 shows 
the variation with temperature of this quantity, for the same parameters 
wpA = 10, A/A = 20 as in Figs. 2-4 (we have taken as our unit scale A the 
value for the odd system (A)=201, which differs from the even system value 
by less than 0.1%). 

In the BCS case, for a finite system in grand-canonical eqnilibrinm, the 
free energy Fbcs{{N)) is a fnnction of the expectation value (A) which may 
vary continuously from (A) = Aq = 200 to (A) = Aq-I-I = 201. On the other 
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Figure 3: Entropy of an odd-number system for = 10, wpA = 200 
and (N) = 201. The solid curve corresponds to the variational solution with 
projection on an odd number of particles. Around T = O.lA it has a Sakur- 
Tetrode form and it tends to ln2 for T —0. The dashed curve (BCS) is 
obtained with the standard hnite-temperature BCS formalism. The dashed- 
dotted and dotted curves are obtained by dropping the interactions. The 
latter unprojected non-interacting case coincides with the entropy of the 
Fermi gas. 
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Figure 4: Specific heat of an odd-number system for wpA = 10, wpA = 200 
and (N) = 201. Symbols and abbreviations are the same as in Fig. 3. The 
solid curve reaches the classical value for one degree of freedom near T = 0.1 A 
while the BCS value tends to zero, as exhibited by the inset. The projected 
and unprojected curves with no interaction are identical and have the linear 
behaviour of the Fermi gas. 


96 





hand, dFBcs/d{N) = fi{{N)) is the BCS chemical potential as fnnction of 
(N) at the considered temperatnre. Accordingly, we can obtain the difference 
SFbcs = -^BCs(-^o + 1) ~ -^BCs(^o) by integration of /r((iV)) between Nq and 
Nq + 1. In the present model with A"o + 1 eqnidistant twofold degenerate 
single-particle levels, for symmetry reasons, /i(Ao -|- 1) is eqnal to the energy 
of the level at mid-spectrnm (level 101) which we take eqnal to zero. The 
Fermi level fi{{N)) lies halfway between two consecntive levels for even (N), 
and coincides with a level for odd (N) provided the occupation number is 
close to 1 (resp.O) at the bottom (resp. top) of the active space. Moreover, 
when wpT S> 1 or wpA S> 1, /i((iV)) is expected to vary smoothly, without 
oscillations. Hence we hnd 

a(W)-^(W-A'o-1) , (6.94) 

from which we derive the BCS difference of free energies: 

rNo+l ^ ^ 1 

SFpcs = / /i((iV))d(iV) ~ . (6.95) 

JNo 4:WF 


At this level of approximation one hnds that 5 Fbcs is temperature inde¬ 
pendent, a result supported by the numerical evaluation, as illustrated by 
Fig. 5. One observes moreover on Fig. 5 that for temperatures above 0.4A 
the projected and BCS curves are identical. 

At lower temperatures, the projected free-energy difference 6F = F-{Nf)+ 
1 ) — Fj_(Ao) is seen to decrease almost linearly from the zero-temperature 
value which, by extrapolation, can be estimated to lie close to A. Let us 
see how this behaviour, found numerically, can be understood from the ap¬ 
proximation (|6.74| ) for the grand potential difference SFq, derived in the 
low-temperature regime (|6.67|) that is relevant here. To obtain Eq. (|6.74|) we 
approximated the solution of the full self-consistent equations by a pertur¬ 
bative treatment starting from the BCS solution. In order to relate 6F to 
6Fq = Fq_{^) — FG_|_(/r), we note that the derivative with respect to /i of the 
approximate expression ( |6.74|) is practically zero. (Indeed dA/dfi\^=o = 0 
for symmetry reasons.) This means that the relation between {N) and fi 
is the same for the even-A and odd-A projected solutions. Hence this re¬ 
lation is also the same as for BCS (Eq. (|6.94|) ). Since F = Fq + fJ,{N) 
and since fi vanishes for (A) = Aq -|- 1, the difference 6Fq coincides with 
F_(Ao + 1) — E+(Ao -|- 1). The extra contribution F^{Nq -I- 1) — F+(Ao) to 
5F is evaluated by integration of dF^/d{N) = p((A)) as in Eq. (|6.95|) . We 
hnd: 

1 IT 

6F = SFq — - -~ A — - -ln( 87 rA tel T) . (6.96) 


4wf 


4wp 
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Figure 5: Difference between the free energies of neighbouring odd {{N) = 
201) and even {{N) = 200) systems as function of temperature for the same 
parameter values wpA = 10 and wpA = 200 as for Figs. 2-4. The solid and 
dashed curves correspond to self-consistent calculations with and without 
projection on number-parity, respectively. Above the crossover temperature 
Tx = 0.3A, the even-odd effect disappears, as for Ap in Fig. 2. Below Tx the 
projected difference 6F decreases nearly linearly. It lies close to the dashed- 
dotted curve which displays the low-temperature perturbative approximation 









This expression is represented on Fig. 5 by the curve labelled “proj.pert.”. It 
is very close to the full solution of the projected equations up to T = 0.3A. 

Thus, our results agree with the observation PT[] that, when the spacing 
between single-particle levels is small compared to the gap {wpA ^ 1), 
number-parity effects disappear above some temperature Tx which lies below 
the critical temperature. Actually, Fig. 5 shows that for wpA = 10 a sharp 
crossover takes place around Tx = 0.3 A (that is, around 0.5 Tc), between the 
A-parity dependent and independent regimes. Below this value, projection 
strongly modifies 6F, while parity effects disappear above. Figs. 2, 3 and 4 
exhibit the identity of the various BCS, even-A and odd-A projected results 
above T = 0.4A. 

More generally, for arbitrary values of wpA much larger than 1, the differ¬ 
ence 5F is given at low temperature by the approximation (|6.96|) , at higher 
temperature by the value —l/{4w-p) associated with the same trivial parity 
effect than for a non-interacting system or for the BCS approximation. The 
crossover between the two regimes takes place around the solution of the 
equation 

^ (6.97) 


A 


IntcpAi 


/SttTx 


This crossover temperature is also visible on Fig. 2, which shows that Af(T) 
is the same for BCS and for both the projected cases in the range < T < 
Tc , where Tx ~ 0.3A for wpA = 10. 

When the size of the system decreases, wpA also decreases. The crossover 
temperature Tx thus rises; from Eq. (|6.97 ) one hnds that it reaches the critical 
temperature Tc = 0.57A for wpA ~ 1.5. Thus, if the distances between levels 
become comparable to the gap, one can expect that odd-even effects will take 
place over the entire range of temperatures where pairing occurs. Physical 
systems which satisfy this condition seem presently available since one can 
prepare and study ultrasmall aluminium grains Q. 


6.5.4 Ultrasmall Metallic Grains 

We have chosen, as an example of modelisation of ultrasmall metallic grains, 
the parameter values wpA ~ 1.1 and tcpA = 100. This is obtained for 
wpG = 0.23, and for an active space having 101 twofold degenerate levels. 
For the odd (resp. even) system we take (A) = 101 (resp. 100). In Fig. 6 we 
present the gap Ap as function of the temperature for odd and even systems. 
The critical temperature associated with the even-A projection comes out 
larger than the BCS value. Even more striking is the behaviour of the odd- 
A projected curve. Then the gap is zero at T = 0 but takes a non-zero 
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value above some first transition temperature up to a second one, lower than 
the BCS critical temperature. This peculiar behaviour can be explained as 
follows. At T = 0, the unpaired particle occupies one of the two single¬ 
particle states Co = /i with a probability exactly equal to one. This blocks 
the formation of pairs utilizing the two degenerate states with energy cq 
and creates an effective single-particle spectrum with a distance between the 
levels at the Fermi surface which is twice as large as A. Consequently, pairing 
is inhibited. When the temperature grows, the bachelor particle statistically 
populates other single-particle states, as indicated in the discussion above on 
the entropy ( p.87|) . This allows pairing correlations on the strategic level at 
the Fermi surface (cq = p), provided the temperature is not too low. One 
thus observes a reentrance effect, with two transition temperatures : one at 
which the pairing switches on and, at a higher temperature, the usual one at 
which it switches off. 

The entropy of the even system (A)=100 is shown in Fig. 7. Projection 
strongly reduces the entropy below the critical point. This projection effect 
is seen to be more important for the interacting system. At small T, the 
four entropies vanish exponentially. As expected the decrease rate of the 
projected entropy is the fastest. For T > 0.3A, the grand-canonical entropy 
of the non-interacting system agrees with the linear form of Eq. ( b.93|) . It 
bends down below T = 0.3A under the effect of the spectrum discreteness; 
in the domain T > 0.3A, projection again induces only a lowering by ln2. 

Fig. 8 displays the entropy of the odd system (A)=101. When the pair¬ 
ing interactions are present, the BCS formalism, which only constrains the 
average number of particles, yields the result 0 instead of the right zero- 
temperature limit ln2. This correct limit is obtained by the odd-A projec¬ 
tion. When there is no interaction, for T > 0.3A, the grand canonical value 
(dotted curve) is again described by Eq. ( |6.93|) . However the entropy tends 
to ln4 at small T. This comes from the fact that, at T = 0, the Fermi energy 
exactly coincides with the energy of the last occupied level whose average oc¬ 
cupation is equal to 1/2. When this value is inserted in the standard formula 
for the grand-canonical entropy of a non-interacting fermion system (while 
taking into account the twofold level degeneracy), one obtains ln4 instead of 
the expected value ln2. Projection on odd-A parity decreases the entropy by 
ln2 and restores the correct limit. Note that this is not a trivial result since 
the exact particle number is not restored by the projection; only its parity 
is. On Fig. 8 the reentrance phenomenon manifests itself by the fact that the 
entropy of the interacting system is larger than that of the non-interacting 
system over the temperature interval 0.14 < T/A < 0.34; otherwise the effect 
of interactions is rather small. 

The original features associated with values of wpA close to one are even 
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Figure 6: Comparison of BCS and number-parity projected gaps at the Fermi 
surface versus temperature for wpA ^ 1.1, wpA = 100, {N) = 100 (upper 
part) or 101 (lower part of the hgure). Such parameters are illustrative of 
extremely small superconducting grains. The critical temperature is raised by 
even-A^ projection, lowered by odd-A^ projection. In the latter case, pairing 
disappears below a new, lower critical temperature 0.15A, thus displaying a 
reentrance phenomenon. 
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Figure 7: Entropy of an even-number system for wpA ^ 1.1, wpA = 100, 
(TV) = 100. The symbols and notations are the same as in Fig. 3 . Projection 
strongly reduces the entropy. 
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Figure 8: Entropy of an odd-number system for (N) = 101 and the same 
parameters wpA ^ 1.1, wpA = 100 as in Figs. 6 and 7. The symbols and 
notations are the same as in Fig. 3. The projection lowers the entropy above 
Tc, but raises it at lower temperatures with the limit ln2 at T —0. 
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more visible on the difference of the specific heats of neighbouring odd and 
even systems. In Fig. 9 one sees that the projected curve shows discontinuities 
for the even-system critical temperature, and also for the two odd-system 
critical temperatures. The comparison with the same curve with wpA = 10 
shown in the inset of Fig. 9, which corresponds to the second derivative of 
the curves on Fig. 5, displays the striking qualitative change arising from the 
variation of wpA over one order of magnitude. 


6.5.5 Comments 


From our numerical analysis, the pairing properties of heavy nuclei (Sect. 6.5.1) 
appear to be intermediate between that of metallic mesoscopic islands (Sect. 6.5.2) 
and that of ultrasmall aluminium grains (Sect. 6.5.4). 

In both the BCS and projected results, the sharp transition from the su¬ 
perconducting to the normal solution arises from the mean-held nature of the 
approximations. The transitions should be smoothed by a projection onto 
a well-dehned particle number, that we have not performed. This projec¬ 
tion would probably have a negligible inhuence on the results obtained for 
mesoscopic metallic islands where {N) is large and hence has small relative 
huctuations. In nuclei, on the other hand, the smoothing of the supercon¬ 
ducting to normal phase transition due to the hnite number of neutrons and 
protons is experimentally observed in the dependence on angular momen¬ 
tum of the moments of inertia. The same phase-transition smoothing effect 
should also be observable in metallic grains when they become sufficiently 
small, through experiments on temperature-dependent properties or on mag¬ 
netic response. The general equations derived in Sect. 5.2 and 5.3 could be 
helpful to perform the required projection. They could also help to account 


for fluctuation effects, such as those invoked in Ref. 120 
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Figure 9: Difference 6Cy between the specific heats of neighbour odd and 
even systems (TV) = 101 and (N) = 100 (multiplied by the ratio of the 
zero-temperature BCS gap A over T) as function of temperature for the 
same parameters wpA ^ 1.1, wpA = 100 as in Fig. 6-8. The solid and 
dashed curves correspond to self-consistent calculations with and without 
projection on particle-number parity, respectively. The critical temperature 
for the even system and the two critical temperatures for the odd system 
show up as discontinuities. The curve differs qualitatively from that of the 
inset, which for (TV) = 200 and (TV) = 201 corresponds to the parameters 
wpA = 10, wpA = 200 of Figs. 2-5. In this inset, the dashed, dashed-dotted 
and dotted curves are indistinghishable from the temperature axis. 
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7 Summary and Perspectives 


A general variational method for evaluating the properties of hnite systems at 
thermodynamic equilibrium has been proposed and applied to hnite fermion 
systems with pairing correlations. Extensions of the standard mean-held 
theories were thus obtained for two types of problems: 1) the variational 
evaluation not only of thermodynamic quantities, but also of expectation 
values, huctuations or correlations between observables and more generally of 
characteristic functions; 2) the restoration of some constraints on the physical 
space of states or/and of some broken symmetries violated by the variational 
Ansatze. 

To reach this double goal, we relied on an action-like functional (p.4|) 
which is designed to yield as its stationary value the relevant characteristic 
function. Two dual trial quantities enter this functional. The operator I?(m) 
is akin to a statistical operator evolving from 1 to the actual density operator 


D = 


(7.1) 


as u goes from 0 to /3, while the operator A{u) is akin to the observable A 
dehned in Eq. ( |1.6|) . Such a doubling of the number of variational degrees 
of freedom is distinctive of a general class of variational principles built to 
optimize some quantity of interest dehned by given constraints. Here the 
constraint is the (symmetrized) Bloch equation (|2.3|) satished by 'D[u), while 
A{u) appears as a Lagrangian multiplier associated with this constraint. For 
unrestricted trial spaces, the stationarity of the functional (|2.4|) expresses 
the Bloch equation (|2.3| ) for T>{u)] in addition, it yields the equation ( |2.10| ) 
for A{u) which is the Bloch analogue of a (backward) Heisenberg equation. 
These two equations are decoupled and duplicate each other. However, when 
the trial spaces are restricted (as is unavoidable for any application), the two 
equations generally become coupled, with mixed boundary conditions involv¬ 
ing both ends of the integration range [0, j3\. The fundamental thermody¬ 
namic relations are preserved by the variational approximations under rather 
general conditions on the choice of the trial spaces (Sect. 2.2); these condi¬ 
tions are satished by the Ansatze that have been used throughout the present 
article. For conserved observables, the characteristic functions, and therefore 
all the cumulants, are obtained from the approximate partition function via 
a simple shift of the operator K (Sect. 2.3). Moreover, in the limit when the 
sources entering the characteristic function vanish, and provided the trial 
spaces satisfy two additional conditions stated in Sect. 2.4, our variational 
principle reduces to the standard maximization principle for thermodynamic 
potentials. 
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Beginning with Sect. 3, we used the functional (|2.4|) to study systems 
of fermions in which pairing correlations are effective. For the variational 
operators F’(m) and A{u) we choose the trial forms 


f'M(M) = exp(-iM(ti) - 
= exp(— 


(7.2) 

(7.3) 


where ot denotes the 2n operators of creation a\ or annihilation ax; the trial 
quantities are the scalars and /[“I and the 2n x 2n matrices and £[“ 1 . 
The generalized Wick theorem then allows to write explicitly the resulting 


functional as Eq. (|3.22|) . The stationarity conditions are expressed by the set 
of coupled equations (|3.37|) , ( p.43D , ( |3.46| ) and (|3.47|) which constitute an ex¬ 
tension of the Hartree-Fock-Bogoliubov approximation. The equation (|3.37|) 
proceeds forwards with respect to the imaginary time u while Eq. (|3.43| ) 
evolves backwards, since the boundary conditions are hxed at m = 0 for the 
former and sX u = (5 for the latter. 

In Sect. 4 thermodynamic quantities, fluctuations and correlations have 
been evaluated via the expansion of the coupled equations in powers of the 
sources The thermodynamic quantities are given by the zeroth order 
term in the expansion while expectation values (Q 7 ) of single-quasi-particle 
observables ^ 

+-OL^ Q^ol (7.4) 

are given by the hrst order. In both cases, one recovers (Sect. 4.1) the 
standard HFB results in agreement with the discussion of Sect. 2. 

The fluctuations and correlations = \{Q'yQ& + Q8Q^) — {Q-y){Q 5 ) are 
given by the second order. In the case where the two observables Q 7 and Qs 
commute with the operator K, we obtain (Eq. (|4.27 )) 


i s,: F„-' : Q, . (7.6) 

The 2n x 2n matrices and Qs dehned in ( [T4|) are regarded in Eq. (fTSi) 
as vectors in the Liouville space, and the colon symbol stands as a scalar 
product in this space, or equivalently as half the trace (tr 2 n ) in the original 
space of the 2n x 2n matrices (Sect. 4.2.2 and Eq. ([A.l| )). In the Liouville 
space, the stability matrix Fq is dehned in Eq. ([4.20|) as the second derivative 
with respect to the reduced density 77 of the HFB grand potential Fq around 
its minimum, the HFB solution TZq. If the observables and Qs do not 
commute with K, Eq. ( |7.5|) expresses their Kubo correlation (Sect. 4.3.1). It 
then appears as a generalization to the non-commutative case of the Ornstein- 
Zernike formula, here derived from a general variational principle. 
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The ordinary correlation of two observables Q-y and Qs that do not com- 
mnte with K is covered by the more general formula (Eq. (|4.63 )) 


= is, : (K„cothi,3K„)F„-‘ : Qs 


(7.6) 


The matrix Kq is the usual RPA kernel written in the Liouville space (Eq. ([4.47|) 
or ( |A.15|) ). The HFB stability matrix Fq is related to the dynamical RPA 
matrix Kq through Eq. (|4.49|) , or ( |A-22|) . 

Notwithstanding the independent-quasi-particle nature of the operators 
'D{u) and .4,('u), we have obtained non-trivial approximate expressions for 
the correlations. In particular, formulas (|7.5|) or ( [f.6D incorporate long-range 
effects. This was made possible because, through the optimization of the 
functional (|2.4|) , the variational quantities V{u) and A{u) acquire a depen- 


7 - 


deuce on the sources 

In spite of the coupling inherent to the stationarity equations (|3.37|) and 
(|3.43| ), we have found the explicit expressions ( [7.5D and (|7.6| ) for the corre¬ 
lations. Moreover, we have found a simple w-dependence for the solutions 
of these equations to first order in the sources. Despite the compactness of 
the final results (|7.5| ) or ( |7.6| ), their derivation was complicated by the ex¬ 
istence of zero eigenvalues in the kernel Fq, due to some broken invariance, 
and in the kernel Kq (Sect. 4.3.2). Indeed, for the formulas (|7.5|) and ([T6|) 
to be meaningful, one must specify how to handle these zero eigenvalues. 
We showed in Sects. 4.2.3 and 4.3.2 that the precise meaning to be given 
to these formulas depends on the commutation, or non-commutation, of the 
observables Q-^ and Qs with the conserved single-particle operator {N in our 
example) associated with the broken invariance. 

Analyzed in perturbation theory, the formulas ( |7.5|) and ([7.6| ) were shown 
to correspond to a summation of bubble diagrams (Sect. 4.3.3). Let us stress 
that the RPA kernel comes out from our variational approach without any 
additional assumption beyond the choice (|7.2|) , (|7.3|) of the trial Ansatze. In 


this way, the RPA acquires a genuine variational status (albeit not in the 
Rayleigh-Ritz sense). 

Sects. 5 and 6 are devoted to the second type of problems posed in the 
Introduction, namely to those situations which require projections. Though 
the operator K, and hence the exact density operator o, commute with 
any conserved quantity, the corresponding invariance can be broken by the 
variational approximation if the commutation is not preserved by the trial 
operator D(m). In this case, a projection P over an eigenspace of the con¬ 
served quantity is required on both sides of T> to restore the symmetry. On 
the other hand, even when the invariance is not broken, a projection is needed 
if the trial operators act in a space which is wider than the physical Hilbert 
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space Til. (For instance we consider below the problem of extending the 
thermal Hartree-Fock approximation to a system with a well dehned particle 
nnmber.) Sects. 5 and 6 cope with both sitnations. 

In the hrst case where the invariance is broken, we replace the trial ex¬ 
pression (|7.2|) for the density operator by 


V{u) = Pf^‘^{u)P 


(7.7) 


where still has the independent quasi-particle form ( [f.2|) ; for A{u) we 

keep the Ansatz m- In order to handle (O) we take advantage of the fact 
that the projection P can be written as a sum of group elements which 


have an independent quasi-particle form, so that (|7.7|) is also a sum of terms 
of the T-type. The resulting functional (Eq. (|5.16|) ) and coupled stationarity 
equations ( (|5.22|) , (^.23|) , (|5.27|) and (|5.28|) ) are written in Sects. 5.2 and 5.3. 
Although these coupled equations are noticeably more complicated than the 
unprojected ones, the number of variational parameters, and therefore the 
dimension of the numerical problem, do not increase. 

The case of unbroken P-invariance, where the trial operators T^^{u) and 

f[a] 

(u) commute with the group elements is discussed in Sect. 5.4. Since 
P now commutes with the Ansatz o reduces to 




(7.8) 


which entails several simplifications in the action-like functional (Eq. ( b.4(j| )) 
and in the stationarity equations. We have examined in some detail (Sect. 5.4.2) 
the limiting case A(,^) = 1 which admits the explicit solution (|5.43|) for P{u) 
and A{u). This leads to further simplihcations in the functional (Eq. ( b.50| )) 
and in the self-consistent equations ( b.51j - 5.54 ). Several thermodynamic con¬ 
sistency properties were verified. As an example the projected entropy, en¬ 
ergy and grand potential were shown to satisfy the usual thermodynamic 
identities (Eqs. ( |2.26| - |2.29|) ). 

In Sect. 6 the theoretical framework elaborated in Sect. 5.4 was applied 
to study the effects of the particle-number parity in a finite superconducting 
system at thermodynamic equilibrium. In the trial density operator ([7.8|) , 

P given by (|5.3|) is the projection on either even or odd particle number. 
This problem is of interest for recent or future experiments on mesoscopic 
superconducting islands, small metallic grains or heavy nuclei. The explicit 
form of the coupled equations (see Eq. (|6.9| )) was written in Sect. 6.1. The 
parity-projected grand partition function (or more generally the character¬ 
istic function for conserved single-particle operators QP was evaluated in 
Sect. 6.2, where Eq. (|6.23|) , together with by the definitions ( |6.25| - b.2^) , was 
shown to replace the standard HFB self-consistent equation. 
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Sect. 6.3 deals with a system (electrons in a snpercondnctor or nncleons 
in a deformed nncleus) governed by a BCS model Hamiltonian, where the 
non-interacting part involves twofold degenerate levels p with energies , 
and the interaction takes place only between the two states of each level. 
The BCS gap is now replaced by the quantities dehned in Eq. (|6.50|) , 
and the usual BCS equation is replaced by the number-parity projected self- 
consistent equation ( |6.56D , 


A 






L^rjq ' U ^q 

^(€, - + AJ 1 + 'll ro 


(7.9) 


in which Gpq is a pairing matrix element of the BCS Hamiltonian ( 5.34 ). The 
notation tp stands for 

fp = tanhi/3ep , (7.10) 

where the quasi-particle energy Cp, given by Eq. (|6.54|) , differs from the stan¬ 
dard BCS expression ^(cp — /i)^ -|- Ap. The number ro is dehned (Eq. ( |6.46| )) 
as 


I'd= n 


(7.11) 


The formulae (|7.9| - |7.1lD depend on the parity of the particle number N 
through r] = -|-1 (even N) or rj = —1 (odd N). The effects of the pro¬ 


jection appear indirectly through the occurence of Cp in (|7.10|) , and directly 
through the last fraction of Eq. O)- As discussed in Sect. 6.3.3, this frac¬ 
tion is larger than one for even-A systems, smaller for odd-A systems. As 
a consequence, the even-A projected gap is larger than the BCS gap, while 
the odd-A projected gap is smaller, in agreement with the idea that pairing 
is inhibited in small odd-A systems. The analysis of this fraction reveals, 
moreover, that the odd-A projection differs from BCS more than the even-A 
projection. 

In Sect. 6.4 we considered some limiting situations. At zero temperature 
the BCS and even-A projected gaps are equal. However, while the BCS 
quantities deviate from their zero-temperature limit by terms of the form 
e“^^, the deviations are of the form for the even-A projection. The 

differences are more striking between BCS and the odd-A projection. At zero 
temperature, the BCS gap is larger than the odd-A projected gap. When the 
temperature grows, the latter starts to increase before reaching a maximum. 
The odd-A entropy tends to ln2 when T tends to zero, reflecting the twofold 
degeneracy of the ground state; in Sect. 6.5.3 we commented upon the fact 
that this result is not completely trivial. At low temperature the odd-A 
entropy coincides with the Sakur-Tetrode entropy of a single-particle with a 
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mass A moving in a one-dimensional box, which is consistent with gapless 
elementary excitations in odd-iV systems. When the single-particle level 
density is sufficiently large, the critical temperatures become equal for the 
BCS and the odd or even projected solutions. Moreover for a large system, 
or when the level spacing becomes small with respect to the pairing gap at 
the Fermi surface, Eq. (|7.9|) reduces to the usual BCS equation. However, 
the value of the projected (odd or even) entropy is lower than the BCS value 
by ln2. 

In order to push further the comparison between the BCS and the pro¬ 
jected results we presented in Sect. 6.5 three schematic calculations which 
exemplify three different physical situations. The crucial parameter turns out 
to be the product wpA where A is the zero-temperature BCS gap and wp the 
single-particle level density at the Fermi surface. This parameter wpA, in¬ 
terpreted as a “number of Cooper pairs”, indicates how many single-particle 
levels lie within the energy range of the gap. 

When wpA is large, the projection effects are weak. Fig. 2 illustrates 
the case wpA = 10. The BCS and even-iV projected gaps are almost undis- 
tinguishable. The relative difference between the BCS and odd-iV gaps is 
l/2tCFA at T = 0. The difference between the free energies of the odd and 
even systems can be estimated perturbatively from the BCS solution. The 
outcome (Eq. (|6.96| )), equal to A — l/Awp at T = 0, decreases almost linearly 
with T and displays a cross-over towards the asymptotic value —l/dtcp at a 
temperature Tx smaller than the critical temperature (see Fig. 5). Between 
Tx and T^, the gap function is the same wathever the A-parity. These re¬ 
sults are consistent with experiments on mesoscopic superconducting islands 
which show that, when wpA S> 1, the odd-even effects disappear at some 


temperature below TdlM 


The value wpA ~ 2 corresponds roughly to the case of a superdeformed 
heavy nucleus (Fig. 1). In the even system, the difference between the pro¬ 
jected and BCS gaps becomes more appreciable than in Fig. 2. In the odd 
system, the projected gap at zero temperature is reduced by about 30% with 
respect to BCS. This result can also be obtained from a BCS calculation in 
which the pair formation is blocked for the level that coincides with the Fermi 
surface. As the temperature grows so does the odd-A projected gap. The 
BCS, even and odd projected gaps merge here again at higher temperature 
where they all display the same critical behaviour. 

Fig. 6, where wpA ~ 1.1, corresponds to a situation which could be 
encountered in ultrasmall Aluminium grains)^. In the even-A projected 
state, the critical temperature is higher than for the BCS state. In the odd 
system a new phenomenon occurs: there is no gap at zero temperature. 
However, the phase diagram (T, A) obtained by projecting on odd number 
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exhibits a reentrance effect : when the temperature increases, a gap appears 
at a hrst transition temperature, reaches a maximum and disappears at a 
second transition temperature smaller than the BCS one. This behaviour can 
be understood as follows : at T = 0, at least one particle remains unpaired 
and it fully occupies one of the two single-particle states at the Fermi surface, 
forbidding the formation of a pair with components on this strategic level. 
When the temperature increases, the unpaired particle partially frees up 
this level, which then becomes available for the formation of Cooper pairs. 
According to Fig. 9, the existence of two transitions in odd ultrasmall metallic 
grains could be experimentally detected by specific heat measurements of odd 
and even systems. 

Besides its relevance for hnite superconducting systems, one can imagine 
using the variational method introduced in this paper for studying other 
physical problems. In particular the formulation of Sect. 5.4, which treats 
situations with unbroken P-invariance, is applicable in various problems. 
Consider for instance a finite fermion system at equilibrium without pairing 
correlations. In this case the exponents in the trial quantities (E3) and 
(|7.3| ) need only contain single-particle terms in a|aj which commute with N. 
The projection P has the form (|5.2|) . A variational method to evaluate the 
canomca/partition function is then furnished by Sect. 5.4.2, which constitutes 
a projected extension (with a well-defined number of particles) of the thermal 
Hartree-Fock approximation. It suffices in Eqs. (|5.43| - ^~44D to take a trial 
operator Hq of the single-particle type. This procedure also provides, along 
the lines of Sect. 2.3, a consistent variational estimate for the characteristic 
function of any conserved single-particle observable (commuting with both N 
and H). For other single-particle observables (commuting with N only) one 
should resort to the method of Sect. 5.4.1 where the trial operators (u) and 
Tkl (u) have a more complicate w-dependence. Likewise, for hnite systems 
without odd-multipole deformations, the projection over a given space-parity 
can be performed by implementing the operator (|5.6|) in Sects. 5.4.1 and 5.4.2. 

In case the symmetry is broken (Sects. 5.2 and 5.3), the functional (|5.16|) 
takes a more intricate form, due to the occurence of two projections. One 
could then make use of a further approximation by using the scheme described 
in Ref. where an integral such as ( ^.21) on the group elements is replaced 
by a truncated sum. This would generate approximate expressions for the 
canonical thermodynamic functions and correlations which should be easier 


to handle numerically than the general expressions (|5.23| - |T^) and ( ^ 

TH). 

In the variational methods described above, one could imagine imple¬ 
menting trial quantities 'P(m) and A{u) with a more general form than (|7.2|) 
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and (|773|). For instance Ref. has demonstrated, in the related context of 
time-dependent Hartree-Fock eqnations, the feasibility of extensions where 

(u) is mnltiplied by some polynomial in creation and annihilation opera¬ 
tors. Using the methods developed here, one conld obtain extensions of the 
HFB formalism, or of the projected HFB formalism, which wonld take into 
acconnt a wider class of correlations. 

The variational setting of Sect. 2 conld easily be adapted to problems 
of classical statistical mechanics. A first possibility would be to start from 
the classical version of the functional (p.4|) in which traces, density operators 
and observables are replaced by integrals, distribution functions and ran¬ 
dom variables in phase space. The commutative nature of these quantities 
implies that a characteristic function would be obtained as a maximum. Al¬ 
ternatively, one could start from the quantal functional (|2.4|) and then apply 
a semi-classical limiting process, for instance through the use of the Wigner 
transform. The latter procedure would yield in addition the hrst quantum 
corrections to the classical limit. This could be of interest for heavy nuclei 
or for sufficiently large clusters of atoms and molecules. 

Minor modifications are needed to deal with systems of interacting bosons. 
Besides the change of the anticommutation rules (|3.2| ) into commutation 
ones, linear terms in a and should be added to the observables (|3.12|) and 
to the exponents of the trial operators (|7.2| ) and ( |7.3|) . This should result 
in a systematic variational approach of Bose-Einstein condensation in hnite 
interacting systems, a phenomenon which has become of direct experimental 
interest with the recent achievement of condensates in dilute atomic gases 
at ultra-cold temperatures. It is worthwhile to note that here again the 
projection on the right number of particles is an important requirement. 

Another physical problem where a projection is essential arises in particle 
physics, when mean-field like approximations violate the colour invariance. 
One may imagine restoring it by means of a projection on colour singlets 
along the lines of Sect. 5. 

In conclusion we would like to stress the consistent character of the present 
method, besides its generality and its flexibility. We have verified the agree¬ 
ment of our approximation scheme with thermodynamics. We have recov¬ 
ered naturally the RPA (the complete series of bubble diagrams). Finally, 
although we have not touched upon the subject in this already long article, 
we point out the complementarity of the present static approach with the 
variational treatment of dynamical problems of Ref. |^. Indeed, the use 
of the Bloch equation as a constraint for the initial state combines coher¬ 
ently with the use of the (backward) Heisenberg equation as a dynamical 
constraint; a variational principle comes out which (within restricted trial 
spaces) optimizes both the initial state and the evolution of the system. 
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This comprehensive static and dynamic variational method should provide a 
convenient tool for evaluating time-dependent correlations, and hence trans¬ 
port properties, either in the grand-canonical formalism as in Sect. 4 or with 
projections as in Sects. 5 and 6. 
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A Geometric Features of the HFB Theory 


A.l The Reduced Liouville Space 


The Liouville representation of quantum mechanics (see for instance 0, 
Sect. 3) relies on the idea that the algebra of observables, when regarded 
as a vector space, can be spanned by a basis of operators. Any observable 
is then written as a linear combination of the operators of the chosen basis, 
with coefficients interpreted as (covariant) coordinates. Cor relatively, for any 
state, the expectation values of the operators of the basis are regarded as the 
(contravariant) coordinates of this state. In the HFB context, it is sufficient 
to consider the subalgebra of quadratic forms of creation and annihilation 
operators a\ (equal to either Oj and a^i). With the notation (|3.12|) , keeping 
aside the trivial c-number q, we are led to regard the coefficients Qw as 
the covariant coordinates of the operator \ ot\Q\\ia\i. The factor | 
and the constraint 




imposed on Q are consistent with the duplication 
a\ax' = -oa'^I + ^aa' = - o-A'Ma^a/.'VA + '^AA' occuring in the operator 

basis. We shall thus consider Qaa' both as a matrix in the 2n-dimensional 
A-space of creation and annihilation operators (Sect. 3.1) and as a vector in 
the 2n x 2n-dimensional reduced Liouville space of observables. 

Likewise, since in the HFB approximation we focus on density operators 
of the form (|3.4|) and since they are characterized by the contraction matrices 
7l\y defined by Eqs. ( p.8| ) and (|3.1CI| - |3.1ID , the contravariant coordinates of 
these density operators reduce to the set T^aa', with the constraint (|3.9|) . 
Again, IZxy is both a matrix in the A-space and a vector in the reduced 
Liouville space of states. 

As in the full space of observables and states, the expectation value of an 
operator in a state characterized by its contractions TZxy, is given 

here by the scalar product of the vectors IZ and Q in the Liouville space. We 
denote this scalar product as 


Q : TZ — ^tT2n QTZ — I X! 2aa''^a' 


(A.l) 


AA' 


In (|A.1|) the colon symbol : indicates therefore both a summation on a 


twisted pair of indices and a multiplication by the factor A accounting for 
the duplication of coordinates both in 7Z and Q. 

Among others, we shall use in this Appendix the following two matrices 
in Liouville space : 

l(AA')(At/i') = 2 (5 a/^'5^a' , S(AA')(/i/i') = 2 cxA/^a^/A' . (A.2) 

They are invariant under the exchange (AA') (/u/u') (symmetry in Liouville 

space). One has also = 1. Moreover, the relation (^.9|) implies that the 
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variation 61Z of a contraction matrix satisfies 


S : (57^ = (T57^V = -57^ 


(A.3) 


A.2 Expansion of the HFB Energy, Entropy and Grand 
Potential 

The HFB energy, defined by E{TZ} = Tr VK/Tr f) (inclnding therefore 
the term —fi{N)), is given by Wick’s theorem as the fnnction ( p.31D of the 
Liouville vector TZw. Its partial derivatives have been dehned in Eqs. ( p.38 - 
3 .401) as the 2n x 2n matrix Tl{TZ} constrained by (|3.41|) . In the Lionville 
space, in agreement with ( fl.38| ), Hw appears as the gradient of E, that is, 
a covariant vector which, with the notation (|A.1|) , reads 


5E{7^} = : 57^ . 

Using (|3.7|) and ( p. 8 |) , the HFB entropy is evalnated as 

S{n} =-Ft V\nV /Ft V - lull V 

= —^tT 2 n [HhilZ + (1 — 7^)ln(l — TV) ] 


(A.4) 


(A.5) 


where the two terms nnder the trace are eqnal. Similarly to (|A-4|) , its gradi¬ 
ent, written by means of (|3.9|) so as to satisfy the same identity (p.41|) as Ft, 
is given by 

6S{n} = C:6n , (A.6) 

with C = [ln(l — Tl)/TZ]. 

The second derivatives of E{71} and dehne the twice covariant 

tensors E and S entering the expansions 

E{n + 6n} = E{n} + H{n}:6n + ^6n:E:6n , (a.7) 

S{n + 671} = 5{77} + C:6n+ \57Z : S : 577 + .. . . (A.8) 

Both matrices E and S are symmetric (for instance, E(aa')(/^^') = E(^^/)(aa'))- 
Moreover they satisfy 


S;E;S = E, S:S:S = S 


(A.9) 


so that E : 67Z obeys the same symmetry relation (|A-3|) as 67Z. The explicit 
form of the second derivative E of E{7Z} comes ont from ( |3.31|) and 577. = 
E : 57Z. When 7Z\\i and 67Zx\i are hermitian matrices in the A-space, the 
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matrix S of the second variations of SjT?.} is always negative in ( |A.8| ). As 

, (A.IO) 


shown in it can be expressed as 

1 




= -2 


dn 


X/i' 


n + v{i-n) 


fiX' 


n + v{i-n) 

which, in a basis where 71 is diagonal (TZxy = Tlx ^aa'); rednces to 

Cx-C^ 


)(AA')(m/^') 


= 2 6x,i'S. 


XiM'Of^X' 


Tlx — Til, 


(A.ll) 


The ratio is eqnal to —1/Tlx{^ — Tlx) when Tl^ = Tlx- 

The standard HFB variational principle (Sects. 4.1) provides the grand 
potential at grand canonical eqnilibrinm as the minimnm of 


Fg{K} = B{K}-1s{R} 


(A.12) 


which is attained for Tl = TIq. From (|A-4|) and (|A-6|) , the stationarity condi¬ 
tion yields the self-consistent eqnation 


«{Ko} = i£i, 


(A.13) 


which dehnes the HFB solution Tlo, in agreement with Eq. (|4.1|) . Around 
this equilibrium, the expansion of the grand potential Fq{TI} reads 


Fq{TIq 6Tiy — Fq^TIq} \5Tl ; Fq : 5Tl -|- ... 
where Fq = Eq — 4 Sq is a positive matrix, like — Sq. 


(A. 14) 


A.3 The HFB Grand Potential and the RPA Equation 

The expansion of Eqs. (|3.52|) around = TIq (or around Tl^^°^ = TIq) leads 
to Eqs. ([4.461) which involve the RPA kernel Kq dehned in Eq. ([4.47]) as 


Ko : = [7^{7^o} , 5TI] + 


dH{Tl} 


/ dllaa' 

flfl' FF 




(A.15) 


'1Z='7Zq 


In Liouville space the matrix Kq is non symmetric ; it depends on Tlo and on 
the parameters of K. We show here that Kq is directly related to the matrix 
Fq which enters (|A.14|) , or (|4.20|) . 

To this aim, we hrst introduce in the Liouville space a new matrix which 
generates the commutators with Tl. Its action on any vector M. is dehned by 


C : M = -M : C = [Tl, M] 


(A.16) 
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Accordingly, its expression as a matrix is 

C(aa')(mm') = , (A.17) 

where the last two terms yield the same contribntion as the hrst two when 
applied to a vector M. satisfying ( |A.3D (i.e., S : Ad = —Ad). The matrix C 
is antisymmetric, and it moreover satisfies 

C : S = S ; C = -C . (A.18) 

In a basis where 77 is diagonal, it rednces to 


^{xx'){^l^l') = 2 (^A/iA^A'(77 a — 77^) , 


(A.19) 


as it is obvions from (|A.16|) 

The commutators with the matrix C can then be represented in Lionville 
space by the prodnct C : S which acts according to 


C : S : 577 = S : C : 577 = -577 : C : S = [£, 577] 


(A.20) 


Indeed, in a basis where 77 is diagonal, this identity, analogous to ( |A.16|) , is 
a direct consequence of (|A.11|) and (|A.19| ). Note that the matrices C and S 
commute. 

The hrst term of ( A.15 ) is thus proportional to ( A.20 ), in which C and 
S have to be taken at the stationary point where the relation £0 = lii(l ~ 
'^o)/'^o = /37d{77o} is satished. To write the second term, we note that 






(A. 21 ) 


a consequence of (|A.7|) . Hence, using (|A.16|) and then adding (|A.20|) , we hnd 
for any 77 and 577 : 

= -C ; E ; 577 + : S : 577 


fi/i 


= -C ; F : 577 


The sought for expression of the kernel ([A.15|) is then simply 

Ko = -Co:Fo , 


(A.22) 

(A.23) 


where the matrices Co and Fq are evaluated at the equilibrium point 77 = TZq. 
The RPA kernel Kq is thus directly related to the matrix Fq 0 / the second 
derivatives of the grand potential. 
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At the minimum (attained at TZo) of the trial grand potential ^^{7^}, the 
matrix Fq is non-negative. Therefore the eigenvalues o/Kq, which are the 
same as those of —: Cq : F^^ in the space of Liouville vectors satisfying 
e), are real. We thus recover the well-known consistency between the 
stability of the variational HFB state and the nature of the associated RPA 
dynamics in real time. Moreover, due to the antisymmetry of Co, the non¬ 
vanishing eigenvalues of Kq come in opposite pairs. 


A.4 Riemannian Structure of the HFB Theory 


By relying on the very existence of the von Neumann entropy S{D), it has 
been shown that the set of density operators D can be endowed with a natural 
Riemannian structure]^. The basic idea is to introduce a metric such that 
the distance ds^ between two neigbouring states D and D -|- 6D is dehned by 


ds^ = -S^S 


(A.24) 


With this metric the general projection method of Nakajima-Zwanzig 
appears as an orthogonal projection in the space of states. A reduction of 
the exact description through projection on a subset of simpler states, in our 
case the HFB states of the form ( |3.4| ), induces from ([A.24| ) a natural metric 
on this subset. It is thus natural to dehne — S, the positive symmetric matrix 
describing the second variations of the entropy ( |A.3| ), as a metric tensor in 
the reduced 2n x 2n Liouville space of states. This dehnition allows us to 
introduce distances between neighbouring HFB states, and also to establish 
a correspondance between the covariant and contravariant components of 
the vectors of Sect. A.l. More precisely, this canonical mapping relates the 
variations of the observables Cat to the variations of those states which 
are their exponentials; it reads : 


S:5TZ = 5C 


(A.25) 


We can also introduce the twice contravariant metric tensor —in the 
subspace of single quasi-particle observables by inverting S in Liouville space 
( : S = 1 ). In a basis where TZ and C are diagonal, its matrix elements 

are 

- 7 ?, _ - 7 ? 

(A.26) 


q-l — ^ ^ 

A 


-'A 

where the ratio is —7 ^a(1 — TZ\) when 7?.^ = 1Z\. In an arbitrary basis, it can 
be written as 


j-i 

’(AA')(m/.') 


= -2 


du 


..uC 


o'- + 1 


A/i' 


g(l-u)£ 

-I- 1 


(A.27) 


/iA' 
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Finally, we can introduce as in (|4.7|) the reduced thermodynamic potential 


lnC{/l} = InTr exp 


— OL^ Col 
2 


= -tr 2 „ln(l + e 


(A.28) 


The expression of the hrst-order term in the expansion in powers of 5C, 


!!<{/: + 5C} = InCj/:} -n:5C- -5C:S-^ :5C + ... 


(A.29) 


exhibits the fact that the entropy ( A.8 ) is the Legendre transform of the re¬ 
duced grand potential (|A-28|) , while the second-order term is seen to describe a 
distance between two neighbouring observables. This last feature is the coun¬ 
terpart in the reduced HFB description of a structure already recognized ||40[ 
in the full Liouville space. 

All these properties, on which we shall not dwell here, confer a geometric 
nature to the reduced HFB description and set it into the general projection 
method 


A.5 Lie-Poisson Structure of the Time-Dependent HFB 
Theory 

We shall now introduce another geometric type of structure, of a symplectic 
type. This will allow us to regard the reduced energy E{71} or the grand 
potential Fq{TZ\ as a classical Hamiltonian for the dynamics associated with 
the HFB theory. 

The evolution of the variables T^aa' is governed in the time-dependent 
mean-field HFB theory by the equation 

ih^ = [n{n }, 7^] , (A.30) 

written in the 2n x 2n matrix representation. 

In order to analyze this equation, we hrst note that the basic operators 
aai of the Liouville representation are the generators of a Lie group with 
the algebra 


[cKAai/, 


OL^OL^, 


11 _ 


t 




(A.31) 


Following the arguments of Ref.[39|, where a similar formulation was set up 


for the Hartree-Fock case, we note that the contractions TZxy are through 
(|3.8| ) in one-to-one correspondence with the operators q;aQ:^/. Regarding the 


TZxxds as a set of classical dynamical variables, we introduce among them a 
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Poisson structure; indeed, the Lie-algebra relation (|A.31|) suggests to charac¬ 
terize the Lie-Poisson bracket between the dynamical variables by the Poisson 
tensor 

ih {TZw , P'/j./i'} = C(AA'){At/i') 5 (A.32) 

where each C is the linear function of the TZ^s inferred from the right hand 
side of ( |A-31|) - As a matrix C is readily identihed with the matrix introduced 
in ( A.16 . A!T7 ) which described, in the Liouville space, the commutation with 
TZ. While the metric of Sect.A.4 was characterized by the symmetric Riemann 
tensor — S, the present Poisson structure is characterized by the antisymmet¬ 
ric Poisson tensor C. For more details on Poisson structures, see for instance 
Ref.g. 

In contrast to the standard Poisson brackets between canonically conju¬ 
gate position and momentum variables, which are c-numbers, the Lie-Poisson 
brackets (|A.32|) depend on the dynamical variables TZw. In the Liouville rep¬ 
resentation, the bracket of two functions / and g of the 7^’s is evaluated by 
saturating their gradients with the tensor C according to 


■h If I 


(A.33) 


The compatibility of (|A-33|) with the dehnition (|A-32| ) results from the prop¬ 
erty (|A.18|) of the tensor C and from the relation 




(A.34) 




which itself is a consequence of (|3.8|) . 

Using the rules ( |A.33|) and ( |A.34|) , we identify the r.h.s. of (|A.30|) as a 
Lie-Poisson bracket since the time-dependent mean-held HFB equation reads 


d7^ 

dt 


that is 


rr|W{K},7J] 

ih dTZ 


1 BE 
2ih dTZ 


:C:(1-S) = 


1 BE . ^ 
ih BTZ BTZ 


^ = {E{K],K} 


(A.35) 


(A.36) 


This equation has thus the form of a classical dynamical equation of motion ; 
the HFB energy E{TZ} plays the role of a classical Hamiltonian, while a Lie- 
Poisson structure is associated with the dynamical variables TZ by Eq. ( |A.32|) . 

Since the gradient (|A-6|) of S'!??.}, as a matrix in the A-space, commutes 
with TZ, the relation (|A.lti|) implies that the Lie-Poisson bracket {S'!??.} , TZ} 


121 



























vanishes. Hence, we can alternatively regard the grand potential (|A.12|) , 
instead of E{TZ}, as a classical Hamiltonian and write 


^ = {Fam , K] 


(A.37) 


This remark is nsefnl when the time-dependent mean-field eqnation (|A.3(]| ) 
is used in the small amplitude limit to describe approximately collective ex¬ 
citations around the grand-canonical equilibrium state. By expanding IZ as 
77o -l- 577, Eq. ( |A.30|) then reduces to the RPA equation associated with the 
HFB approximation, 

^ = 4ko:577 , (A.38) 

dr in 

where the RPA kernel Kq is expressed by (|A.15|) . We can now identify the 
r.h.s. of (|A-38|) with the first-order term {6Fq{TZ} , 77o} in the expansion 
of (|A.37| ) in powers of 577 around TZq. Using the expansion (|A.14|) and the 
expression (|A-32|) of the Lie-Poisson bracket, we find that (|A-38|) is equivalent 
to 

7hd577/dt =-Co : Fo : 577 , (A.39) 

and thus recover the expression ( |A.23| ) for Kq. This new derivation of ( |A.23|) 
provides us with a simple interpretation of the RPA equation (|A-39|) as an 
equation for small motions in ordinary classical dynamics. Indeed, small 
changes of canonically conjugate variables qi,Pi around a minimum of a 
Hamiltonian H{q,p} are governed by 




{q, q} 

{p, q} 


{7, P} 

{p, p} 


d^H d'^H 


d^H 



dpdq dpdp 

(A.40) 

where one recognizes the same multiplicative structure as in the r.h.s. of 
(|A-39|) , within the replacement of the ordinary Poisson bracket by (|A-32|) 
and of H{q,p} by Fg{77}. 

For small deviations around an arbitrary TDHFB trajectory, the time- 
dependent RPA kernel K is defined by 


zh ^ = K : 577 ^ [nm , 577] + ^^577^^., 77] 

at , UrCnnl 


(A.41) 


Since we no longer lie at the minimum of Tg{77}, this variation (|A.41|) in¬ 
volves not only the expansion of Fq but also of 77, which appears both 
directly and through the Poisson tensor C. We thus get, in the product 


K : 577 = Pf : 5C - C : E : 577 

= -C ; F : 577 + (R: - p-^C) : 5C 


(A.42) 
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an additional term proportional to the deviation between [3T-i{7l} and C = 
ln[(l — TV)/TZ\. Apart from this term, which arises from the dependence of 
the strncture constants C on the variables TZ, the present classical dynam¬ 
ics differs from an ordinary Hamiltonian dynamics through the presence of 
vanishing eigenvalues in the tensor C. The Poisson structure generated by 
(|A.32|) therefore differs from the usual symplectic structure associated with 
canonically conjugate variables because some combinations of the dynami¬ 
cal variables TZ have a vanishing bracket with any variable. This property 
reflects the existence of structural constants of the motion, which never vary 
whatever the effective Hamiltonian in ( |A.3fi|) . In particular, the eigenvalues 
of 7Z and S{7Z} always remain constant. 
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B Liouville Formulation of the Projected Finite- 
Temperature HFB Equations 


Using the Liouville-space formalism introduced in Appendix A, it is possible 
to recast the coupled equations (^.231) and (|5.28|) governing the evolution of 
and dehned in Eqs. ( |3.6|) and (|3.16 - |3.17|) , so as to make these 

equations formally similar to Eqs. (p.37|) and (|3.43|) . 

To this aim, we hrst introduce two operators and whose action 
on a Liouville vector Q is given by 


■ Q — 7"H ^ 

rpM] . Q 

= 7"kl QT'M ^ 

(B.l) 

They satisfy the relations 




S : T^^^ : I] = T^“^ , 

S ; T^'^^ 

; S = T^'^' , 

(B.2) 

rp[a]“l _ rp[a]'^ 

pk]“^ — 

p[d]'^ 

(B.3) 


T' 

where the symbol ^ denotes the matrix transposition in the Liouville space. 
Next, we dehne the matrix 


M 


[.d.a] _ 

(AA')(/./.') = 


+ {n^ad9'a] _ jll.d.a]^^^^^j^ldg'ag] _ ^ 


and the two Liouville vectors 

g[.d.a] ^ J ^^gg'_ 'j^[gdg'a]'j'j^^j^lgdg'a]jj^lgdg'a] 

-|_(£'{7^M9'“]} _ ^(da)jf^J^lgdg'a] _ -^[.d.a]^j 
g[a.d.] ^ j ^^gg'j^(^'j^[agdg']'j^^'j^[agdg']jf-^_'j^[agdg']'j 

_|_(£'{7^k9<^3']} _ ^(adjjf^y^lagdg'] _ 'j^[a.d.]'jj 


(B.4) 


(B.5) 


After multiplication to the right by T[“l, Eq. (|5.23|) can be reformulated 


as 


rlT'H 1 ^ _ 

]y[[-'^-«] ■ (_ _T'M a _|_ (gl-d-a] _|_ rp[a] ^ . gla.d.]\ _ g 

du 2^ 

By introducing in the Liouville space the matrix 

j^[a.d.] ^ rp[a] . j^[.d.a] . rp[d] 

namely 


(B.6) 


(B.7) 


(B.8) 
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and the vectors 


-^[.d.a] _ jyj-[.d.a]-l . g[.d.a] 

we can rewrite 


n 


a.d.] 


= ]y[[a.rf.]-l . g 


a.d.] 


as 


a form which resembles closely ( |3.37|) . 

In the same way Eq. (|5.28|) can be transformed into 

HTbl 1 

_ _ — _|_ q-[a]j_^[d.a.]-j 


(B.9) 


(B.IO) 


(B.ll) 


where the matrices “ and “ l in the single-particle space are obtained 
from Eqs. (1130) and (|B.8| - |B.9|) by the exchanges a <-> d and g <-*• g'. 
One major difference, however, between the conpled eqnations (|B.10| - pn|) 
and Eqs. ( |3.37| ) and (|3.43|) is that the former involve fonr different qnantities 

instead of two for the 

latter. 

The matrices and reqnired for the calcnlation of the qnan¬ 
tities and appearing in ( |B.11| ) are related to and 

throngh 

^[.a.d] ^ ]r^[a.d.]T ^ ^ _ (B.12) 

From the dehnition (|B.4| - |B.5|) and (|B.7D of the matrices and Mb '^-1 

and of the Lionville vectors and ^[“■‘^■1, one establishes the relations 


S ; ; S = , S : : S = 

5] . g[-d.a] _ _ g[.d.a\ ^ . g[a.d.] _ _ gla.d.] 


(B.13) 

(B.14) 


One dednces that, as the rednced HFB Hamiltonian T-[{TZ}, the vectors in Li¬ 
onville space and verify the relation (|3.41|) which in the Lionville 

formnlation becomes : 


5 ] . _ _ J-^^l-d.a] ^ . 'j^la.d.] _ _ 


a.d.] 


(B,15) 


Throngh the exchanges a ^ d and g g' one extends the results (|B.13| - pT5|) 
to the matrices and vectors 

When the operator A (dehned in Eq. (|1.6|) ) is hermitian, Eqs. (|B.10|) , 
( |B.11D , (^.23| ) and (|5.28| ), have hermitian solutions. Indeed, in such a case 
these equations preserve the hermiticity of and and the reality of 
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/[“I and To prove it, one first notices that the relations ^1, 

^ ^[d]t ^ j-[d] 


-^[gdg'ajt _ -^[ag' ^dg p 


(B.16) 


'99 


and (|5.2CI|) implies 


^99'* 

= Z^' ' 

(B.17) 

g- 1 , one 

sees from Eq. (|B.17|^ that Z 

is real. The 

with the 

equality E{Tl^ = E{7lY in 

Eqs. (15.19; 


^[.d.a]t _ -^[a.d.] ^ 

(B.18) 

^{da) 

* _ ^{ad) 

(B.19) 


Finally, on formulae ([B.4| - |B.5[) and (|B.8|) , one readily checks the relations 



[.d.a] * T.y|-[a.d.] 

(A'A)(m'm) “ ’ 

(B.20) 

g[d-a.]i _ g[.a.d] 

gl.d.a]] ^ g[a.d.] _ 

(B.21) 

They yield the equations 



-^[d.a.p _ -^[.a.d] ^ 


(B.22) 

which entail the hermiticity of dT^'^l/dw and dT^^l/dw. 
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